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Abstract: 

This review addresses recent developments in nonequilibrium statistical physics. Focusing 
on phase transitions from fluctuating phases into absorbing states, the universality class 
of directed percolation is investigated in detail. The survey gives a general introduction 
to various lattice models of directed percolation and studies their scaling properties, field- 
theoretic aspects, numerical techniques, as well as possible experimental realizations. In 
addition, several examples of absorbing-state transitions which do not belong to the di- 
rected percolation universality class will be discussed. As a closely related technique, we 
investigate the concept of damage spreading. It is shown that this technique is ambiguous 
to some extent, making it impossible to define chaotic and regular phases in stochastic 
nonequilibrium systems. Finally, we discuss various classes of dep inning transitions in 
models for interface growth which are related to phase transitions into absorbing states. 
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1 Introduction 



Random behavior is a common feature of complex physical systems. Although systems in 
nature generally evolve according to well-known physical laws, it is in most cases impossible 
to describe them by means of ab initio methods since details of the microscopic dynamics 
are not fully known. Instead, it is often a good approximation to assume that the individual 
degrees of freedom behave randomly according to certain probabilistic rules. For this 
reason methods of statistical mechanics become essential in order to study the physical 
properties of complex systems. In this approach a physical system is described by a 
reduced set of dynamical variables while the remaining degrees of freedom are considered 
as an effective noise with a certain postulated distribution. The actual origin of the noise, 
which may be related to chaotic motion, thermal interactions or even quantum-mechanical 
fluctuations, is usually ignored. Thus, statistical mechanics deals with stochastic models 
of systems that are much more complicated in reality. 

A complete description of a stochastic model is provided by the probability distri- 
bution -Pt(s) to find the system at time t in a certain configuration s. For systems at 
thermal equilibrium this probability distribution is given by the stationary Gibbs ensem- 
ble Peq{s) ~ Q-'H{s)/kBT ^ where 'H(s) denotes the microscopic Hamiltonian Q. In principle, 
the Gibbs ensemble allows us to compute the expectation value of any time-independent 
observable by summing over all accessible configurations of the system. However, in most 
cases it is very difficult to perform the configurational sum. In fact, although numerous 
exact solutions have been found Q, the vast majority of stochastic models cannot yet be 
solved exactly. In order to investigate such nonintegrable systems, powerful approxima- 
tion techniques such as series expansions and renormalization group methods Q have 
been developed. Thus, in equilibrium statistical mechanics, we have a well-established 
theoretical framework at our disposal. 

From the physical point of view it is particularly interesting to investigate stochas- 
tic systems in which the microscopic degrees of freedom behave collectively over large 
scales Collective behavior of this kind is usually observed when the system undergoes 

a continuous phase transition. The best known example is the order-disorder transition in 
the two-dimensional Ising model, where the typical size of ordered domains diverges when 
the critical temperature is approached |^. In most cases the emerging long-range cor- 
relations are fully specified by the symmetry properties of the model under consideration 
and do not depend on details of the microscopic interactions. This allows phase transi- 
tions to be categorized into different universality classes. The notion of universality was 
originally introduced by experimentalists in order to describe the observation that several 
apparently unrelated physical systems may be characterized by the same type of singular 
behavior near the transition. Since then universality became a paradigm of the theory of 
equilibrium critical phenomena. As the number of possible universality classes seems to 
be limited, it would be an important theoretical task to provide a complete classification 
scheme, similar to the periodic table of elements. The most remarkable breakthrough 
in this respect was the application of conformal field theory to equilibrium critical phe- 
nomena ||8-1C], leading to a classification scheme of continuous phase transitions in two 
dimensions. 
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In nature, however, thermal equilibrium is rather an exception than a rule. In most 
cases the temporal evolution starts out from an initial state which is far away from equi- 
librium. The relaxation of such a system towards its stationary state depends on the 
specific dynamical properties and cannot be described within the framework of equilib- 
rium statistical mechanics. Instead it is necessary to deal with a probabilistic model for 
the microscopic dynamics of the system. Assuming certain transition probabilities, the 
time-dependent probability distribution Pt{s) has to be derived from a differential equa- 
tion, the so-called Fokker-Planck or master equation. Nonequilibrium phenomena are also 
encountered if an external current runs through the system, keeping it away from thermal 



equilibrium |11|. A simple example of such a driven system is a resistor in an electric 
circuit. Although the resistor eventually reaches a stationary state, its probability dis- 
tribution will no longer be given by the Gibbs ensemble. As a physical consequence the 
thermal noise produced by the resistor is no longer characterized by a Gaussian distri- 
bution. Similar nonequilibrium phenomena are observed in catalytic reactions, surface 
growth, and many other phenomena with a flow of energy or particles through the system. 
Since nonequilibrium systems do not require detailed balance, they exhibit a potentially 
richer behavior than equilibrium systems. However, as their probability distribution can- 
not be expressed solely in terms of an energy functional 'H(s), the master equation has to 
be solved, being usually a much more difficult task. Therefore, compared to equilibrium 
statistical mechanics, the theoretical understanding of nonequilibrium processes is still at 
its beginning. 

The simplest nonequilibrium situation is encountered if a single or several particles in 
a potential are subjected to a random force. Important examples are the Kramers and 
the Smoluchowski equations describing the evolution of the probability distribution for 
Brownian motion of classical particles in an external field (for a review see jjl2| ) . But even 
more complicated systems, for example one-dimensional tight-binding fermion systems as 
well as electrical lines of random conductances or capacitances, can be described in terms 
of discrete single-particle equations [|^. A more complex situation, on which we will focus 
in the present work, emerges in stochastic lattice models with many interacting degrees 



of freedom. A well-known example is the Glauber model [14| which describes the spin 



relaxation of an Ising system towards the stationary state. The corresponding master 



equation in one dimension was solved exactly by Felderhof [15|, who mapped the time 
evolution operator onto a quantum spin chain Hamiltonian that can be treated by similar 
methods as in Ref. llf 



One motivation for today's interest in particle hopping models originates in the study 



of superionic conductors in the 70's |17,1^. In the superionic conductor Agl, for example, 
the Ag"*" ions may be viewed as particles moving stochastically through a lattice of I~ 
ions. Each lattice site can be occupied by at most one Ag"^ ion, i.e., the particles obey an 
exclusion principle. It was observed experimentally that the conductivity of Agl changes 
abruptly when the temperature is increased, indicating an underlying order-disorder phase 
transition of Ag+ ions. Assuming short-range interactions, this phase transition was ex- 
plained in terms of a model for diffusing particles on a lattice Subsequently, particle 
hopping models have been generalized to so-called reaction-diffusion models by including 
various types of particle reactions or external driving forces It should be noted that 
particles in a reaction-diffusion model do not always represent physical particles. More- 
over, the reactions are not always of chemical nature. For example, in models for traffic 



flow individual cars are considered as interacting particles piy. Similarly, electronic ex^ 
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citations of certain polymer chains may be viewed as particles subjected to a stochastic 
temporal evolution p2| . 

The dynamic properties of a reaction-diffusion model on a lattice are fully specified 
by its master equation p^-p5|. In a few cases it is possible to solve the master equa- 
tion exactly. During the last decade there has been an enormous progress in the field 
of exactly solvable nonequilibrium processes. This development was mainly triggered by 
the observation that the Liouville operator of certain (l-l-l)-dimensional reaction-diffusion 
models is related to Hamiltonians of previously known quantum spin systems. For exam- 



ple, as first realized by Alexander and Holstein |26|, the symmetric exclusion process can 
be mapped exactly onto the Schrodinger equation of a Heisenberg ferromagnet. This type 
of mapping was extended to various other one-dimensional reaction-diffusion processes by 



Alcaraz et al. |27|, allowing exact methods of many-body quantum mechanics such as the 



Bethe ansatz and free-fermion techniques to be applied in nonequilibrium physics |25 
Moreover, novel algebraic techniques have been developed in which the stationary state 
of certain reaction-diffusion models is expressed in terms of products of non-commuting 
algebraic objects [pst| . 

In spite of this remarkable progress, the majority of reaction-diffusion models cannot 
be solved exactly. It is therefore necessary to use approximation techniques in order to 
describe their essential properties. The oldest approximation method is the law of mass ac- 
tion, where the reaction rate of two reactants is assumed to be proportional to the product 
of their concentrations. This mean field approach is justified if diffusive mixing of particles 
is much stronger than the influence of correlations produced by the reactions. Mean- field 
techniques have been applied successfully to a large variety of reaction-diffusion systems. 
The study of pattern formation in nonlinear reaction-diffusion models, for example, is 
essentially based on a mean- field approach |3^. However, as has already been realized 
by Smoluchowski fluctuations may be extremely important in low-dimensional sys- 
tems where the diffusive mixing is not strong enough For example, if particles of 
one species diffuse and annihilate by the reaction A -\- A — > 0, the standard mean- field 
approximation predicts an asymptotic decay of the particle concentration as p{t) ^ t^^. 
In one dimension, however, the density is found to decay as p{t) ~ This slow decay 
is due to fluctuations produced by the dynamics, leading to spatial anticorrelations of the 
particles. The existence of such fluctuation effects has been confirmed experimentally by 
measuring the luminescence of annihilating excitons on polymer chains 



As in equilibrium statistical mechanics, nonequilibrium phenomena are particularly 
interesting if the system undergoes a phase transition, leading to a collective behavior of the 
particles over long distances. There is a large variety of phenomenological nonequilibrium 
phase transitions in nature, ranging from morphological transitions of growing surfaces IS^] 



to traffic jams [34|. It turns out that the concept of universality, which has been very 
successful in the field of equilibrium critical phenomena, can be applied to nonequilibrium 
phase transitions as well. However, the universality classes of nonequilibrium critical 
phenomena are expected to be even more diverse as they are governed by various symmetry 
properties of the evolution dynamics. On the other hand, the experimental evidence for 
universality of nonequilibrium phase transitions is still very poor, calling for intensified 
experimental efforts. 

In the present work we will focus on nonequilibrium phase transitions in models with 
so-called absorbing states, i.e., configurations that can be reached by the dynamics but 
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cannot be left. The most important universality class of absorbing-state transitions is 
directed percolation (DP) [^]. This type of transition occurs, for example, in models for 
the spreading of an infectious disease. In these models the lattice sites are considered as 
individuals which can be healthy or infected. Infected individuals may either recover by 
themselves or infect their nearest neighbors. Depending on the infection rate, the spreading 
process may either survive or evolve into a passive state where the infection is completely 
eliminated. In the limit of large system sizes the two regimes of survival and extinction are 
separated by a continuous phase transition. As in equilibrium statistical mechanics, the 
critical behavior close to the transition is characterized by diverging correlation lengths 
associated with certain critical exponents. Similar spreading processes with the same 
exponents can be observed in models for catalytic reactions, percolation in porous media, 
and even in certain hadronic interactions. It turns out that all these phase transitions 
belong generically to a single universality class, irrespective of microscopic details of their 
dynamic rules. In view of its robustness, the DP class may therefore be as important 
as the Ising universality class in equilibrium statistical mechanics. Amazingly, directed 
percolation is one of very few critical phenomena which cannot be solved exactly in one 
spatial dimension. Although DP is easy to define, its critical behavior is highly nontrivial. 
This is probably one of the reasons why DP continues to fascinate theoretical physicists. 

The present review addresses several aspects of nonequilibrium phase transition^. In 
the following Section we introduce elementary concepts of nonequilibrium statistical me- 
chanics such as the master equation, reaction diffusion processes, Monte Carlo simulations, 
as well as the most important analytical methods and approximation techniques. The 
third Section discusses the problem of directed percolation, including a comprehensive in- 
troduction to DP lattice models, basic scaling concepts, approximation techniques, as well 
as field-theoretic methods. In view of the robustness of DP, it is particularly interesting 
to search for non-DP phase transitions which usually emerge in presence of additional 
symmetries. These exceptional universality classes, which have attracted considerable at- 
tention during the last few years, will be reviewed in Sec. 4. Sec. 5 discusses a simulation 
technique, called damage spreading, which has been used in the past to search for chaotic 
behavior in random processes. It is shown that this technique suffers of severe conceptual 
problems, making it impossible to define chaotic phases. We also discuss the critical be- 
havior of damage spreading transitions which are closely related to phase transitions into 
absorbing states. Finally, we turn to depinning transitions in models of growing interfaces 
which are related to nonequilibrium phase transitions into absorbing states as well. As 
it is not intended to cover the whole field of nonequilibrium critical phenomena, we will 
not address various related topics such as self-organized critical phenomena |36|, modified 
reaction-diffusion processes [32|, the dynamics of reacting fronts |37|, driven diffusive sys- 
tems 1 11], and recent results on spontaneous symmetry breaking and phase separation in 
one-dimensional systems ||3^. For further reading we will give references to related fields. 
Supplementary information concerning the definition of tensor products, the derivation of 
the effective action of Reggeon field theory, Wilson's shell integration, and the one-loop 
integrals for DP are given in the appendices For easy reference we also append a list 
of frequently used symbols and abbreviations. 



^The review is based on a Habilitation thesis submitted by the author to the Free University of Berhn 
in June 1999. 
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2 Stochastic many-particle systems 



In this Section we discuss elementary concepts of nonequilibrium statistical mechanics. In 

order to introduce basic notions, we first consider the example of a simple random walk. 
Turning to many-particle systems we introduce the asymmetric exclusion process which 
is a model for biased diffusion of many particles on a onc-dimcnsional line. Moreover, 
we explain the standard mean field approach to reaction-diffusion processes. In order to 
demonstrate the importance of fluctuations, two simple lattice models with particle reac- 
tions will be discTissed, namely coagulation 2A A and pair annihilation 2A 0. It 
turns out that in one dimension the temporal evolution of these systems differs signifi- 
cantly from the mean field prediction, proving that fluctuations may play an important 
role. Furthermore, we review basic concepts of numerical simulation techniques comparing 
different update schemes. Finally we turn to certain analytical methods by which reaction 
diffusion models can be solved exactly. In particular we discuss a recently introduced 
algebraic technique which allows the stationary state of certain nonequilibrium models to 
be expressed in terms of products of noncommutative operators. 



2.1 The one- dimensional random walk 

In order to introduce basic concepts of nonequilibrium statistical physics, let us first con- 
sider a simple symmetric random walk on a one-dimensional line. The 'configuration' of 
this dynamical system at time t is characterized by the position of the walker s{t). A 
random walk may be defined either on a continuous manifold or on a lattice. If both 
position s and time t are discrete variables, an unbiased random walk may be realized by 
the random process 

sit + l) = sit) + X{t) . (1) 

In this expression X{t) = ±1 is a fiuctuating random variable with correlations 

{X{t)) = , {X{t)X{t')) = 5t,t' , (2) 

where (...) denotes the average over many realizations of randomness. While the individ- 
ual space-time trajectory of a random walker is not predictable, the probability distribution 
Pt{s) to find the walker after t time steps at position s evolves deterministically according 
to the so-called master equation 

P^^^{s) = \{Pt{s-l) + Pt{s + l)) . (3) 

Assuming the particle to be initially located at the origin Po{s) = dg^o-, this difference 
equation is solved by 

R(,) = /^W^)/2) if ^ + ^ even, 

\0 ff t + s odd. 

If space and time are continuous, the motion of a random walker may be described by a 
stochastic Langevin equation 

dts{t) = at) , (5) 
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where, according to the central hmit theorem, is a Gaussian white noise with zero 
mean and correlations (C(*)C(*')) = ^S{t — t'). The Langevin equation may be regarded 
as a continuum version of Eq. (^). Starting from the origin s(0) = the mean square 
displacement of the random walker grows as = Jq dti Jg dt2(C(^i)C(^2)) = ^t. The 

resulting probability distribution 



is a solution of the Fokker-Planck equation ||l^ 



dtPt{s) = ^d',Pt{s) (7) 

which can be seen as a variant of the master equation in a continuum (^). 

The example of a random walk is particularly simple as it involves only one degree 
of freedom. In order to describe systems with many particles, it would seem natural 
to introduce several degrees of freedom si,S2, • • • , where Sn denotes the position of the 
n-th particle. However, this approach is restricted to systems with a conserved number 
of particles. For systems with non-conserved particle number it is more convenient to 
introduce local degrees of freedom for the number of particles located at certain positions 
in space. 



2.2 The master equation 



Stochastic systems with many particles are usually defined on a d-dimensional Euclidean 
manifold representing the physical 'space'. Attached to this manifold are local degrees 
of freedom characterizing the configuration of the system. Depending on whether the 
spatial manifold is continuous or discrete, the local degrees of freedom are introduced 
as continuous fields or local variables residing at the lattice sites. Furthermore, a time 
coordinate t is introduced which may be interpreted as an additional dimension of the 
system. Therefore, stochastic models are said to be defined in d+1 dimensions. Since t may 
be continuous or discontinuous, we have to distinguish between models with asynchronous 
and synchronous dynamics. 



Asynchronous dynamics 

Stochastic models with continuous time evolve by asynchronous dynamics, i.e., transitions 
from a state s into another state s' occur spontaneously at a given rate Wg—ts' > per 
unit time. It can be shown that in the limit of very large systems sizes the temporal 
evolution of the probability distribution Pt{s) evolves deterministically according to a 
master equation with appropriate initial conditions [p3-25|. The master equation is a 



linear partial differential equation describing the flow of probability into and away from a 
configuration s: 

^Pt(s) = Ws'^sPtis') - J2 Ws^s'Ptis) . (8) 



gain loss 

The gain and loss terms balance one another so that the normalization ^^gPtis) = 1 is 
conserved. Since the temporal change of Pt{s) is fully determined by the actual probability 
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distribution at time t, the master equation describes a Markov process, i.e., it has no 
intrinsic memory. Moreover, it is important to note that the coefficients Wg^s' ^ire rates 
rather than probabihties. Thus, they may be larger than 1 and can be rescaled by changing 
the time scale. 

Using a vector notation (see Appendix ^) the master equation (^) may be written as 

dt\Pt) = -C\Pt) , (9) 

where \Pt) denotes a vector whose components are the probabilities Pt{s)- The Liouville 
operator C generates the temporal evolution and is defined through the matrix elements 

{s'\C\s) = -Ws^s' + Ss,s' ^ Ws^s'^ . (10) 

s" 

A formal solution of the master equation is given by \Pt) = exp(—Ct)\PQ) , where |Po) 
denotes the initial probability distribution. Therefore, in order to determine \Pt), the 
Liouville operator has to be diagonalized which is usually a nontrivial task. 

Apart from very few exceptions, stochastic processes are irreversible and therefore not 
invariant under time reversal. Hence the Liouville operator C is generally non-hermitean. 
Moreover, it may have complex conjugate eigenvalues, indicating oscillatory behavior. 
Oscillating modes are not only a mathematical artifact, but can be observed experimentally 
in certain chemical reactions such as the Belousov-Zhabotinski reaction Due to the 
positivity of rates, the real part of all eigenvalues is nonnegative, i.e., the amplitude of 
excited eigenmodes decays exponentially in time. The spectrum of the Liouville operator 
includes at least one zero mode C\Ps) = 0, representing the stationary state of the system. 
Moreover, probability conservation can be expressed as {l\Pt) = 1, where (1| denotes the 
sum vector over all configurations (cf. Appendix ^) . Consequently the Liouville operator 
obeys the equation {1\C = 0, i.e., the sum over each column of £ vanishes. 

Synchronous dynamics 

If the time variable t is a discrete quantity, the model evolves by synchronous dynamics, 
i.e., all lattice sites are simultaneously updated according to certain transition probabilities 
Ps^s' S [0, 1]. The corresponding master equation is a linear recurrence relation 

Pt+l{s) = Pt{s) + Y.Ps'--sPt{s') - Y,Ps^s'Pt{s) , (11) 



gain loss 

which can be written in a compact form as a linear map 

|Pt+i) = T\Pt) , (12) 

where T is the so-called transfer matrix. A formal solution is given by \Pt) = T^\Pq). As 
can be verified easily, the conservation of probability {l\Pt) = 1 implies that (1|T = (1|, 
i.e., the sum over each column of the transfer matrix is equal to 1. 

There has been a long debate which of the two update schemes is the more 'realistic' 
one. For many researchers models with uncorrelated spontaneous updates appear to be 
more 'natural' than models with synchronous dynamics where all particles move simul- 
taneously according to an artificial clock cycle. On the other hand, many computational 
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Fig. 1: The partially asymmetric exclusion process in one dimension with synchronous (sublattice-parallel) 
and asynchronous (random-sequential) dynamics. 



physicists prefer stochastic cellular automata with synchronous dynamics since they can be 
implemented efficiently on parallel computers. However, as a matter of fact, in both cases 
the dynamic rules are simplified descriptions of a much more complex physical process. 
Therefore, it would be misleading to consider one of the two variants as being more 'nat- 
ural' than the other. Instead, the choice of the dynamic procedure should depend on the 
specific physical system under consideration. Very often both variants display essentially 
the same physical properties. In some cases, however, they lead to different results. For 
example, models for traffic flow with synchronous updates turn out to be more realistic 
than random sequential ones. Another example is polynuclear growth (see Sec. |6.3| ), where 
a roughening transition occurs only when synchronous updates are used. 



2.3 Diffusion of many particles: The asymmetric exclusion process 

One of the simplest stochastic many-particle models is the partially asymmetric exclusion 
process on a one-dimensional chain with N sites In this model hard-core particles 
move randomly to the right (left) at rate q ((?~^). An exclusion principle is imposed, i.e., 
each lattice site may be occupied by at most one particle. Therefore, attempted moves are 
rejected if the target site is already occupied. In the following we assume closed boundary 
conditions, i.e., particles cannot leave or enter the system at the ends of the chain. The 
configuration s = {si,S2,... ,S7v} of the system is given in terms of local variables Sj, 
indicating presence (sj = 1) or absence (s, = 0) of a particle at site i. 

Exclusion process with asynchronous dynamics 

Let us first consider the exclusion process with asynchronous dynamics. In this case a 
pair of sites i and i + 1 \s randomly selected. If only one of the two sites is occupied, the 
particle moves with probability q/ {q + q~^) to the right and with probability q~^ /{q + q'"^) 
to the left, as shown in Fig. [l|. Each update attempt corresponds to a time increment of 
l/{N{q + q~^)). Thus, the transition rates are defined by 

Af-l i-l N 

ws^s' = E (n '^^..^0 ( n '^^.--;) ( '^'^^-i ^^«+i'0 ^4+1.1 + , , 

i=l j=l j=i+2 (13) 
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The corresponding Liouville operator can be written as 

7V-1 N~l 

C= ^ 1 1 (g) . . . . . . 1 =: ^ A , (14) 

j-th position 

where 1 denotes a 2 x 2 unit matrix and £j is a 4 x 4 matrix generating particle hopping 
between sites i and i + 1 . In the standard basis (see Appendix ^ this matrix is given by 



/o 
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Vo 








0/ 



(15) 



As can be verified easily, a stationary state of the system is given (up to normalization) 
by the tensor product 

IP.) ^ (;)«(',>... «(v>g(v;). a«) 

Since the vector \Ps) can be written as a tensor product, the local variables Sj are com- 
pletely uncorrelated. Such a state is said to have a product measure. 

As the total number of particles is conserved in the asymmetric exclusion process, the 
dynamics decomposes into independent sectors. In fact, the Liouville operator commutes 
with the particle number operator 

M = Y,m, m=[l I) ■ (17) 



Obviously the vector \Ps) is a superposition of solutions belonging to diff'erent sectors, 
i.e., it represents a whole ensemble of stationary states. To obtain a physically mean- 
ingful solution for a given number of particles, the vector IP,) has to be projected onto 
the corresponding sector. In this sector the stationary system evolves through certain 
configurations with specific weights given by the normalized components of the projected 
vector. 



Exclusion process with synchronous dynamics 

The asymmetric exclusion process with synchronous updates may be realized by introduc- 
ing two half time steps. In the first half time step the odd sublattice is updated whereas 
the even sublattice is updated in the second half time step (see Fig. 0). Note that the use 
of sublattice-parallel updates admits local dynamic rules^. Assuming the number of sites 
N to be odd, the corresponding transfer matrix reads 



where 



T = (T2 r4 % 



N-2) , 



(q + q-^ 



+ 



V 



q 












q + q-^J 



(18) 



(19) 



^In the exclusion process with fuUy parallel updates local moves may overlap, leading to subtle long- 
range correlations (see Ref. pc|]). 
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diffusion coagulation decoagulation annihilation 

Fig. 2: Nearest-neighbor processes in reaction-diffusion models with only one particle species. 

is the local hopping matrix. Again the product state (|l^ is a stationary eigenvector of 
the transfer matrix. Thus, both the asynchronous and the synchronous exclusion process 
have exactly the same stationary properties. 

Asymmetric diffusion in a continuum 

Let us finally turn to asymmetric diffusion on a continuous manifold. In principle it 
would be possible to trace trajectories of individual particles. However, it is much more 
convenient to characterize the state of the system by a density field p{x,t), rendering the 
coarse-grained density of particles at position x at time t. The Langevin equation of such 
a system may be written as 

5<p(x, t) = D[p{x, t)] + U[p{x, t)] + C(x, t) , (20) 

where Z) is a sum of linear differential operators describing spatial couplings, U a potential 
for on-site particle interactions, and (^(x, t) a noise term taking the stochastic nature of 
Brownian motion into account. Since particles do not interact in the present case, the 
potential U vanishes. Moreover, as will be shown below, the noise C(x, t) is irrelevant on 
large scales and can be neglected. Thus, the resulting Langevin equation reads 

9iP(x, t) = ii^s2p(x, t) + i:^9,p(x, t) . (21) 

The second term describes the bias of the diffusive motion which may be eliminated in a 
co-moving frame. Notice that this equation is linear and does not incorporate the exclusion 
principle. In a co-moving frame it reduces to the ordinary diffusion equation. 

The diffusion equation provides a simple example of dynamic scaling invariance. As 
can be verified easily, the equation dtp{'x., t) = D\/'^p{'K, t) is invariant under rescaling of 
space and time 

x^Ax, t^AH, (22) 

where z is the so-called dynamic exponent. Since space and time are different in nature, 
the exponent z is usually larger than 1. The value z = 2 indicates diffusive behavior. 

2.4 Reaction-diffusion processes 

Reaction-diffusion processes are stochastic models for chemical reactions in which parti- 
cles are predominantly transported by thermal diffusion. Usually a chemical reaction in 
a solvent or on a catalytic surface consists of a complex sequence of intermediate steps. 
In reaction-diffusion models these intermediate steps are ignored and the reaction chain 
is replaced by simplified probabilistic transition rules. The involved atoms and molecules 
are interpreted as particles of several species, represented by capital letters A, B,C, . . . . 
These particles neither carry a mass nor an internal momentum, instead the configuration 
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of a reaction-diffusion model is completely specified by the position of the particles. On a 
lattice with exclusion principle such a configuration can be expressed in terms of local vari- 
ables Sj = 0, 1, 2, 3, . . . , representing a vacancy and particles A, B,C, . . . , respectively. 
Sometimes it is even not necessary to keep track of all substances involved in a chemical 
reaction. For example, if a molecule of a gas phase is adsorbed at a catalytic surface, 
this process may be effectively described by spontaneous particle creation — > A without 
modeling the explicit dynamics in the gas phase. Therefore, the number of particles in 
reaction-diffusion models is generally not conserved. 

Apart from spontaneous particle creation — > A many other reactions are possible. 
Unary reactions are spontaneous transitions of individual particles, the most important 
examples being 

A ^ self-destruction, 

A 2A offspring production or decoagulation, 

A ^ B transmutation, 

A ^ A + B induced creation of particles. 

On the other hand, binary reactions require two particles to meet at the same place (or 
at neighboring sites). Here the most important examples are 

2A ^ pair annihilation, 

2 A A coagulation (coalescence), 

A + B ^ two-species annihilation, 

A + B ^ B induced adsorption of particles. 

In addition, the particles may diffuse with certain rates in the same way as in the previ- 
ously discussed exclusion process. A process is called diffusion-limited if diffusion becomes 
dominant in the long-time limit, i.e., the diffusive moves become much more frequent than 
reactions. This happens, for example, in reaction-diffusion models with binary reactions 
when the particle density is very low. On the other hand, if particle reactions become 
dominant after very long time, the process is called reaction-limited. 

In the following we will focus on simple reaction-diffusion models with only one type 
of particles (see Fig. ^). They can be considered as two-state models since each site can 
either be occupied by a particle {A) or be empty (0). Examples include the so-called 
coagulation model in which particles diffuse at rate £), coagulate at rate A and decoagulate 
at rate k: 

A0^0A, AA^A0,0A, A0,0A^AA. (23) 

DA ft 



Using the same notation as in Eq. (|15|), the corresponding nearest-neighbor transition 
matrix is given by 



^coag 



1^ 













D + K 


-D 


-A 





-D 


D + K 


-A 


Vo 


— K 


— K 


2A / 



(24) 



The special property of this model lies in the fact that the empty state cannot be reached 
by the dynamic processes. An exact solution of the coagulation model will be discussed 
in Sec. I2I. 
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Another important example is the annihilation model where particles diffuse and an- 
nihilate: 



A0 ^ 0A 

D 



AA^m 



(25) 



The corresponding interaction matrix is given by 



CI 



annh 











—a\ 





D 


-D 








-D 


D 





Vo 








a / 



(26) 



In the annihilation model the number of particles is conserved modulo 2. As will be shown 



in Sec. 2.8, both models are equivalent and can be related by a similarity transformation. 



2.5 Mean field approximation 

In many cases the macroscopic properties of a reaction-diffusion process can be predicted 
by solving the corresponding mean field theory. In chemistry the simplest mean field ap- 
proximation is known as the 'law of mass action': For a given temperature the rate of 
a reaction is assumed to be proportional to the product of concentrations of the react- 
ing substances. This approach assumes that the particles are homogeneously distributed. 
It therefore ignores any spatial correlations as well as instabilities with respect to inho- 
mogeneous perturbations. Thus, the homogeneous mean field approximation is expected 
to hold on scales where diffusive mixing is strong enough to wipe out spatial structures. 
Especially in higher dimensions, where diffusive mixing is more efficient, the mean field ap- 
proximation provides a good description. It becomes exact in infinitely many dimensions, 
where all particles can be considered as being neighbored. 

The mean field equations can be constructed directly by translating the reaction scheme 
into a differential equation for gain and loss of the particle density p{t). For example. 



in the mean field approximation of the coagulation model (23) the process A 2 A 
takes place with a frequency proportional to Kp{t), leading to an increase of the particle 

density. Similarly, the coagulation process 2A A A decreases the number of particles with 
a frequency proportional to Xp^{t). Ignoring diffusion, the resulting mean field equation 
reads 

dtp{t) = Kp{t) - Xp\t) , (27) 

where k and A are the rates for decoagulation and coagulation, respectively. In contrast 
to the master equation this differential equation is nonlinear. For k > it has two fixed 
points, namely an unstable fixed point at /3 = and a stable fixed point at p = k/\. The 
physical meaning of the two fixed points is easy to understand. The empty system remains 
empty, but as soon as we perturb the system by adding a few particles, it quickly evolves 
towards a stationary active state with a certain average concentration p > 0. This active 
state is then stable against perturbations (such as adding or removing particles). 



The mean field equation (27) can also be used to predict dynamic properties of the 



system. Starting from a fully occupied lattice p{Q) = 1 the time-dependent solution is 
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given by 



In the limit of a vanishing decoagulation rate k — > 0, the two fixed points merge into a 
marginal one. As in many physical systems this leads to a much slower dynamics. In fact, 
for K = Eq. ( p8[ ) turns into 

i.e., the particle density decays asymptotically according to a power law as 

p{t) ~ t-' . (30) 

The stability of the mean field solution with respect to inhomogeneous perturbations may 
be studied by adding a term for diffusion 

9tp(x, t) = Kp{x, t) - Ap2(x, t) + Z)VV(x, t) . (31) 

In the present case the diffusive term suppresses perturbations with short wavelength and 
therefore stabilizes the homogeneous solutions. However, in certain chemical reactions with 
several particle species such a diffusive term may have a destabilizing infiuence. The study 
of mean-field instabilities is the starting point for the theory of pattern formation which 



has become an important field of statistical physics |30|. A very interesting application is 



the Belousov-Zhabotinski reaction |^9| that produces rotating spirals in a Petri dish. 

It may be surprising that even simple reaction-diffusion processes are described by 
nonlinear mean field rate equations, whereas the corresponding master equation is always 
linear. However, mean field and Langevin equations are always defined in terms of coarse- 
grained particle densities involving many local degrees of freedom. These coarse-grained 
densities, which can be thought of as observables in configuration space, may evolve ac- 
cording to a nonlinear laws. A similar paradox occurs in quantum physics: Although the 
Schrodinger equation is strictly linear, most observables evolve in a highly nonlinear way. 



2.6 The influence of fluctuations 

Although the mean field equation (^Tj) includes a term for particle diffusion, it still ignores 
fiuctuation effects and spatial correlations. However, especially in low-dimensional sys- 
tems, fluctuations may play an important role and are able to entirely change the physical 
properties of a reaction-diffusion process. 

In order to demonstrate the influence of fluctuations, let us consider the coagulation 
process 2A — > A. The full Langevin equation for this process reads 

5tp(x, t) = -Ap2(x, t) + L>vV(x, t) + C(x, t) , (32) 

where C(x, t) is a noise term which accounts for the fluctuations of the particle density at 
position X at time t. Clearly, the noise amplitude depends on the magnitude of the density 
fleld p{x,t). In particular, without any particles present, there will be no fluctuations. 
According to the central limit theorem, the noise is expected to be Gaussian with a 
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squared amplitude proportional to the frequency of events leading to a change of the 
particle number. Since the particle number only fluctuates when two particles coagulate, 
this frequency should be proportional to p^(x, i). Following these naive arguments, the 
noise correlations should be given by 

(C(x,t)) = 0, (33) 
(C(x,t)C(x',0) = rp'{^,t)5''{^-^')5{t-t'), 

where T denotes the noise amplitude and d the spatial dimension. The next question 
is to what extent the macroscopic behavior of the system will be affected by the noise. 
Typically there are three possible answers: 

1. The noise is irrelevant on large scales so that the macroscopic behavior is correctly 
described by the mean field solution. 

2. The noise is relevant on large scales, leading to a macroscopic behavior that is dif- 
ferent from the mean-field prediction. 

3. The noise is marginal, producing (typically logarithmic) deviations from the mean- 
field solution. 

In order to find out whether the noise is relevant on large scale we need to introduce the 
concept of renormalization The term 'renormalization' refers to various theoretical 
methods investigating the scaling behavior of physical systems under coarse-graining of 
space and time. Roughly speaking, it describes how the parameters of a system have to be 
adjusted under coarse-graining of lengths scales without changing its physical properties. 
A fixed point of the renormalization flow is then associated with certain universal scaling 
laws of the system. The simplest renormalization group (RG) scheme ignores the influence 
of fluctuations. This approach is referred to as 'mean field renormalization'. Approaching 
the fixed point, the noise amplitude may diverge, vanish or stay finite, corresponding to 
the classification given above. Hence, by studying mean field renormalization, we can 
predict whether fluctuations are relevant or not. 

In the mean field approximation the Langevin equation ( |32[ ) may be renormalized by 
a scaling transformation 

x^Ax, t^AH, p(x,t) ^ A^/)(Ax,A^t), (34) 

where z denotes the dynamic exponent. The exponent x describes the scaling properties of 
the density field itself. If the particles were distributed homogeneously, the field would scale 
as an ordinary density, that is, with the exponent x = ~d- However, in the coagulation 
process nontrivial correlations between particles lead to a different scaling dimension of 
the particle distribution. In fact, invariance of Eqs. (|3^)-(|33|) under rescaling implies that 
z = 2 and x = ~2. Therefore, the noise amplitude scales as 

r ^ A^-'^/^T , (35) 

where d is the spatial dimension. Hence in one spatial dimension fluctuations are relevant 
whereas they are marginal in two and irrelevant in d > 2 dimensions. The value of d where 
the noise becomes marginal is denoted as the upper critical dimension dc- For the coagula- 
tion model the upper critical dimension is dc = 2. Above the critical dimension the mean 
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Fig. 3: Monte Carlo simulation of the (l+l)-dimensional coagulation-difltusion model A + A ^ A with 
random-sequential updates. The figure shows an individual run starting with a fully occupied lattice of 
200 sites. As can be seen, the coagulation process is diffusion-limited. 



field approximation provides a correct description, whereas for d < dc fluctuation effects 
have to be taken into account. This can be done by using improved mean field approaches, 
exact solutions, as well as field-theoretic renormalization group techniques 



A systematic field-theoretic analysis of the coagulation process 2A — > A leads to an 
unexpected result: The noise amplitude F in Eq. ( p3[ ) turns out to be neg ative |4|. 



Consequently, the noise t) is imaginary. This result is rather counterintuitive as we 
expect the noise to describe density fluctuations which, by definition, are real. However, 
since the noise amplitude is a measure of annihilation events, it is subjected to correlations 
that are produced by the annihilation process itself. In one dimension these correlations 
are negative, i.e., particles avoid each other. This simple example demonstrates that it can 
be dangerous to set up a Langevin equation by considering the mean field equation and 
adding a physically reasonable noise field. Instead it is necessary to derive the Langevin 



equation directly from the microscopic dynamics, as explained in Ref. |44| 



2.7 Numerical simulations 

To verify analytical results, it is often helpful to perform Monte Carlo simulations. In order 
to demonstrate this numerical technique, let us again consider the coagulation process 
2A — > A on a one-dimensional chain. For simplicity we assume the rates for diffusion and 
coagulation to be equal. This ensures that particles can move at constant rate irrespective 
of the state of the target site. If the target site is empty, it will be occupied by the moving 
particle. On the other hand, if the target site is already occupied, the two particles will 
coagulate into a single one. Such a move from site i to site j may be realized by the pseudo 
code instruction 

Moved, j) { if (s[i]==l) { s[i]=0; s[j]=l; }}; 

where s [i] denotes the occupation variable Si = 0, 1 at site i. In one dimension particles 
move randomly to the left and to the right. Thus a local update at sites (i, i + 1) may be 
realized by the instruction 
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Fig. 4: Monte Carlo simulation of the coagulation-diffusion model. Left: Decay of the particle density 
p{t) measured in a system with 10^ sites in one, two, and three spatial dimensions, averaged over 10* 
independent runs. The dashed lines indicate the slopes —1/2 and —1, respectively. Right: tp{t) versus t, 
illustrating logarithmic corrections in d = 2 dimensions. 



Update(i) { if (rnd(0 , 1) <0 . 5) Move(i,i+l); 
else Move(i+l,i); }; 

where rnd(0, 1) returns a real random number from a flat distribution between and 1. 
Since the coagulation model evolves by asynchronous dynamics it uses so-called random- 
sequential updates, i.e., the update attempts take place at randomly selected pairs of sites. 
A Monte Carlo sweep consists of N such update attempts: 

for (i=l; i<=N; i++) Udpate(rndint(l,N-l)) ; 

where N denotes the lattice size and rndint (1 , N- 1) returns an integer random number 
between 1 and A^ — 1. Since on average each site is updated once, such a sweep corresponds 
to a unit time step. It can be proven that the statistical ensemble of space-time trajec- 
tories generated by random-sequential updates converges to the solution of the master 
equation (^) in the limit N ^ oo. The above update algorithm can easily be generalized 
to more complicated reaction schemes and higher dimensional lattices. 

The coagulation process with synchronous updates may be simulated by using parallel 
updates on alternating sublattices: 

for (i=l; i<=N-l; i+=2) Update (i); 
for (i=2; i<=N-2; i+=2) Update (i); 

In Monte Carlo simulations most of the CPU time is consumed for generating random 
numbers. Therefore, models with parallel updates are usually more efficient since it is 
not necessary to determine random positions for the updates. In addition, models with 
parallel updates can be implemented easily on computers with parallel architecture ||4^ . 

Fig. ^ shows the particle density as a function of time for the coagulation model 
with random-sequential updates and closed boundary conditions in various dimensions. 
The particle concentration is averaged over 10^ independent runs and plotted in a double- 
logarithmic representation, where straight lines indicate power-law behavior. As expected, 
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the mean-field prediction p(t) ~ 1/t is reproduced in d > 2 dimensions. In one dimension, 
however, the graph suggests the density to decay as 



p{t) ~ (d = 1) 



(36) 



Thus, the simulation result demonstrates that fluctuation effects can change the asymp- 
totic behavior (an exact solution will be discussed in Sec. 2^). At the critical dimension 
d = dc = 2 the density p{t) deviates slightly from the mean-field prediction, indicating 
logarithmic corrections. In fact, as can be shown by a field-theoretic analysis |44], the 
density decays asymptotically as 



Pit) 



for d < 2 
for d = dc 
for d > 2 



(37) 



2.8 Exact results 

Equivalence of annihilation and coagulation processes 

Sometimes it is possible to relate different stochastic processes by an exact similarity 
transformation p^-p8[. For example, the coagulation process 2A A and the annihi- 



lation process 2A — > defined in Sec. 2.4 are fully equivalent if their rates are tuned 



appropriately. More precisely, for a particular choice of the rates it is possible to find a 
similarity transformation U such that 

where we assume the chains to have closed ends, i.e., 

N-l 

^coag ^ ^ ^coag ^annh ^ ^ ^annh (39) 

i=l i=l 

Since £™^s and C'^'^^^ are non-hermitean operators, the similarity transformation U is 
not orthogonal. However, if U exists, the two operators will have the same spectrum 
of eigenvalues. As can be verified easily, the local operators £^°'^^ and Cf^^^ have the 
eigenvalues {0, 0, 2D + k, 2A + k} and {0, 0, 2D, a}, respectively. Therefore, choosing the 
rates 

D = X = l, a = 2, K = 0, (40) 

both operators obtain the same spectrum {0, 0, 1, 1}. Moreover, it can be shown that they 
both obey the same commutation relations, namely the so-called Hecke algebra l49|,^ 

— r. 

CiCiJ^iCi — Ci^iCiCi^i = 2{Ci — £j+i) , (41) 

=0 for >2) . 

Since this algebra generates the spectrum of £, we can conclude that the spectra of 

and £^'1°'^ coincide for an arbitrary number of sites. Obviously, the equivalence of the 

spectra is a necessary condition for the existence of a similarity transformation between 
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Fig. 5: The partially asymmetric exclusion process with particle adsorption and desorption at the bound- 
aries. Left: Particles adsorb at the left boundary at rate a, perform a biased random walk in the bulk 
of the chain, and finally desorb from site at rate /3. Right: Phase diagram for the totally asymmetric 
exclusion process, denoting the low (LD), high (HD) and maximal (MAX) density phases. Analogous 
phase diagrams for partial asymmetry were obtained in Ref. [p7|. 



the two systems. In the present case it is even possible to compute the similarity transfor- 
mation explicitly. It turns out that lA can be expressed in terms of local tensor products 
(see Appendix ^) 



U = u®u® ...®u = i^u = u®^^ , (42) 

i=l 

where 

Consequently the n-point density correlation functions of both models are related by 

{sj^Sj2 • • • Sj^) ^ = 2 l^Sj^Sj2 ■ ■ ■ Sj^) . (44) 
In particular, the particle densities in both models differ by a factor of 2: 

^coag(^^ = 2p^°'^^(t). 

It should be noted that only a subset of initial conditions in the coagulation model can 
be mapped onto physically meaningful initial conditions in the annihilation model (see 
Ref. 0). 

Exact mapping between equilibrium and nonequilibrium systems 

A remarkable progress has been achieved by realizing that certain nonequilibrium models 
can be mapped onto well-studied integrable equilibrium models. More specifically, it has 
been shown that the Liouville operator £ of a nonequilibrium models may be related to the 
Hamiltonian Ti of an integrable quantum spin systems by similarity transformation [l5[| . 
This allows the nonequilibrium model to be solved by exact techniques of equilibrium 
statistical mechanics such as free-fermion diagonalization, the Bethe ansatz, or other al- 
gebraic method s |[2 |. For example, the exclusion process can be mapped onto the XXZ 



quantum chain 1 18, 26, 5^, whereas the coagulation-decoagulation process is related to the 



chain in a magnetic field [53-55|. Exact mappings were also found for higher spin 



analogs |56|. A complete summary of the known results can be found in Ref. [28|. 



21 



In order to demonstrate this technique, let us consider the partially asymmetric exclu- 
sion process in one spatial dimension (see Eq. ^). As will be shown in the following, this 
model is related to the quantum spin Hamiltonian of the ferromagnetic XXZ Heisenberg 
quantum chain with open boundary conditions 7i 



Yji=i where 



Tii 



1 



+ 



q + q~ 



1) + 



q-q 



a: 



(45) 



This quantum chain Hamiltonian generates translations in the corresponding two-dimen- 
sional XXZ model in a strongly anisotropic scaling limit. As can be verified easily, the 
Hamiltonian is non-hermitean for q ^ \. Using the standard basis of Pauli matrices 
(T^ = (5J), 0"^ = (i^*)) and = ( i*! ) the interaction matrix is given by 



q-^ 
-1 

Vo 














0/ 



(46) 



The XXZ Heisenberg chain is integrable by means of Bethe ansatz methods . The 
integrability is closely related to two different algebraic structures. On the one hand, the 
Hamiltonian ( pSj ) commutes with the generators of the quantum algebra Uq[SU{2)] 

(see Ref. |8|) 



KS^K-^ = qS 



± 



K'-K- 
q - q-^ 



where 



N 



K={q 



(47) 



(48) 



k=l 



Roughly speaking, the quantum group symmetry determines the degeneracies of the eigen- 
values of 7i. On the other hand, the generators Tii are a representation of the Temperley- 



Lieb algebra [^| 



-{q + q-^)'Hi , 



m 

-Hi 

[n^,nj] =Ofor \i-j\>2. 



(49) 



This algebra determines the actual numerical value of the energy levels. As realized by 



Alcaraz and Rittenberg [56|, the same commutation relations are satisfied by the transi- 
tion matrix £j of the asymmetric exclusion process in Eq. (|15|). In fact, it can be shown 
that the XXZ chain and the exclusion process with a = /? = are related by a simi- 
larity transformation C = U7iU~^, where U can be written as a tensor product of local 
transformations 



1 
q 
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q^ 



1 

q^ 



N 



1=1 



1 

q' 



(50) 



In order to illustrate how symmetries of the equilibrium model translate into physical 
properties of the stochastic process, let us consider the quantum group symmetry of the 



22 



gain: [0 0000000" » ]0 loss: O O O O O O ^e) » 



[o o o o o o Q cri) » 

Fig. 6: Interparticle distribution functions. The figure illustrates gain and loss processes for the empty- 
interval probability Ie(t) by diffusion (left) and coagulation (right). The same processes take place at the 
left boundary of the interval. 

XXZ model. Regarding the exclusion process this symmetry emerges as a conservation 
of the total number of particles n. The generators act as ladder operators between 
different sectors with a fixed number of particles. The diagonal operator K is proportional 
to q~'^, weighting the sectors as in a grand-canonical ensemble, where q plays the role of 
a fugacity. The partially asymmetric exclusion process is therefore a physical realization 
of a quantum group symmetry with a real- valued deformation parameter q. 

A similar mapping relates the coagulation model 2 A — > A and the ferromagnetic XY 
quantum chain in a magnetic field, which is defined by the interaction matrix 

-Hi = -\ + afaf_,i + af + af^, - 2) . (51) 

In fact, it can be easily verified that the operators TLi satisfy the Hecke algebra (^). The 
chain is exactly solvable in terms of free fermions jl^. The integr ability of the model 
is closely related to a SUq{l\l) quantum group symmetry. In the XFchain this symmetry 
shows up as a fermionic zero mode, leading to two-fold degenerate energy levels. In the 
coagulation model this symmetry emerges as a state without particles that can neither be 
reached nor left by the dynamics. In the (equivalent) annihilation model the symmetry 
appears as a parity conservation law. 

It should be noted that reaction-diffusion models are usually related to ferromagnetic 
quantum chains, the reason being that the diffusion process always corresponds to a fer- 
romagnetic interaction in the quantum spin model. 

Interparticle distribution functions 

Even if a stochastic model can be mapped onto a known equilibrium system by a similarity 
transformation, it is often technically difficult to derive physical quantities such as density 
profiles and correlation functions (6^ . For models with an underlying fermionic symmetry 
an alternative approach has been developed which does not explicitly use a similarity 
transformation. Instead it expresses the state of a model in terms of so-called interparticle 
distribution functions (IPDF) pl|-|63|. Consider, for example, the coagulation-diffusion 
process with asynchronous dynamics on an infinite chain where particles coagulate (^4 + 
A ^ A) and diffuse at unit rates. Let Ii be the probability that an arbitrarily chosen 
interval of £ sites contains no particles. In terms of these empty-interval probabilities 
the master equation can be written in a particularly simple form. Since Ii — I^+i is the 
probability to find a particle at a neighboring site next to the interval of length £, it is 
possible to rewrite diffusion and coagulation in terms of gain and loss processes (see Fig. ^) 

dth{t) = le.iit) - h{t) - h{t) + /,+i(t) (52) 
' r ' ' ' 

gain loss 
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with Io{t) = 1. It is important to note that this particularly simple form requires the 
rates for diffusion and coagulation to be identical. This ensures that the gain processes do 
not depend on whether the target site is already occupied by a particle. If the two rates 
are different, higher-order probabilities for several adjacent intervals have to be included, 
resulting in a coupled hierarchy of equations. The IPDF method exploits the fact that 
this complicated hierarchy of equations decouples for a particular choice of the rates. 

By solving the above equation we can compute the particle density p{t) which is given 
by the probability for an empty interval of length 1 to be absent, i.e., 

p{t) = l-h{t). (53) 

In order to determine the asymptotic behavior of p{t) let us consider the continuum limit 
of Eq. (H) 

{dt-d!)I{l,t)=0 , 1(0,0 = 1, (54) 

where p{t) = deI{i,t)\^^Q. This equation has the solution I{i,t) = 1 — erf(x/\/t). In the 
long time limit, the particle density therefore decays algebraically as 

Pit) ~ , (55) 



confirming the numerical result of Sec. |2.7| . It is interesting to compare this result with the 
mean field approximation (^) which lead to the incorrect result p{t) ~ t~^. Therefore, 
the above exact solution demonstrates that fluctuations may influence the entire temporal 
evolution of a stochastic process. 



The IPDF technique [£l,63| was extended to the coagulation-decoagulation model by 
including the inverse reaction A A + A. Other exact solutions revealed phenomena 
such as anomalous kinetics, critical ordering, nonequilibrium dynamic phase transitions, 
as well as the existence of Fisher waves ||6^-|68|. The IPDF technique was also used to 
study the finite-size scaling behavior of coagulation processes [54, ^|. Even anisotropic 



systems [55| and models with homogeneous or localized particle input |65, 7^] have been 
solved. Nevertheless the IPDF method is a rather special technique which seems to be 
restricted to models with an underlying fermionic symmetry. 



2.9 Experimental verification of fluctuation efl'ects 

The preceding exact calculation proves that the particle concentration in a one-dimensional 
coagulation process decays as p{t) ~ This result differs significantly from the mean- 

field prediction p{t) ~ 1/t. Therefore, the coagulation model provides one of the simplest 
examples where fluctuation effects change the entire temporal behavior of a reaction- 
diffusion process. 

It is quite remarkable that this result could be verified experimentally by analyzing 
the kinetics of laser-induced excitons on tetramethylammonium manganese trichloride 
(TMMC) p2|, [7l|. TMMC is a crystal consisting of parallel manganese chloride chains. 
Laser-induced electronic excitations of the Mn^"'' ions, so-called excitons, migrate along 
the chain and may be interpreted as quasi-particles. The chains are separated by large 
tetramethylammonium ions so that the exchange of excitons between different chains is 
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Fig. 7: Decay of the luminescence of TMMC after excitation by a laser pulse for various energies (figure 
reprinted from Ref. WW)- The straight lines are best fits to Eq. (p9). 



suppressed by a factor of 10^. Therefore, the polymer chains can be considered as one- 
dimensional systems. Because of exciton-phonon induced lattice distortions the motion 
of excitons is diffusive. Moreover, when two excitons meet at the same lattice site, the 
Mn^"*" ion is excited to twice the excitation energy. Subsequently, the ion relaxes back 
to a simply exited state by the emission of phonons. Thus, the fusion of excitons can be 
viewed as a coagulation process 2A A+heat on a one-dimensional lattice. 

The concentration of quasiparticles can be measured indirectly by detecting the lumi- 
nescence intensity I{t) which is proportional to the number of excitons. Eq. ( |52|) predicts 
that I{t) ~ /(0)/(l + at^), where a fixes the time scale and 6 = 1/2. This equation can 
be rewritten as 

= 6logt + loga . (56) 

The experimental results are shown in Fig. |^. The best fits according to Eq. (|56| ) yield 
estimates of about 6 = 0.48(3), being in perfect agreement with the theoretical prediction 
6 = 1/2. In other experiments the polymers chains are confined to small pores. Here 
the excitons perform both annihilation and coagulation processes. The estimates 6 = 
0.55(4) ||72[ and S = 0.47(3) |Q are again in agreement with the theoretical result. 

To summarize, the experimental investigation of excitons on polymer chains confirms 
that the concept of stochastic reaction-diffusion processes is well justified in order to quan- 
titatively predict the behavior of certain complex systems. In addition, these experiments 
prove that fluctuations effects do exist in nature and may change the physical properties 
of the system in agreement with the theoretical prediction. 
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2.10 Dynamic processes approaching thermal equihbrium 

Stochastic dynamic processes also play an important role in the context of equilibrium 
models. As outlined in the Introduction, equilibrium statistical mechanics deals with 
many-particle systems in contact with a thermal reservoir (heat bath) of temperature T. 
In the long-time limit such a system approaches a statistically stationary state where it 
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evolves through certain configurations according to a well-defined probability distribution 
Peq{s). The key property of equilibrium models is the existence of an energy functional Ti. 
associating each configuration s with a certain energy 7i{s). The equilibrium distribution 
Peq{s) is then given by the canonical ensemble |l| 

Peg(5) = |exp(-W(s)/fcBr) , (57) 

where T is the temperature, ks the Boltzmann constant, and Z the partition sum. This 
probability distribution can be used to determine averages of certain macroscopic observ- 
ables by summing over all accessible states. It is important to note that equilibrium 
statistical mechanics does not involve any dynamical aspect. In other words, it is irrele- 
vant how the system evolves through different configurations, one is only interested in the 
relative frequency of certain configurations to be visited in the stationary state. 

Although there is no 'time' in equilibrium statistical physics, one may use dynamic 
random processes as a tool to generate the equilibrium ensemble Peq{s) of a particular 
equilibrium model. More precisely, such a dynamic process evolves into a stationary state 
-foo('S) := limt^oo -ft('S) that coincides with the equilibrium ensemble Peq{s). Generally 
there is a large variety of dynamic random processes that can be used to generate the 
stationary ensemble of a particular equilibrium model. Let us, for example, consider the 
ferromagnetic Ising model on a d-dimensional square lattice 0. Its energy functional is 
given by 

n{a) = -jY, ^i^j > (58) 

<i,3> 

where the sum runs over pairs of adjacent sites, J is a coupling constant, and fij = ±1 
denotes the local spin at site i. The equilibrium ensemble of the Ising model may be 
generated by a dynamic process with synchronous dynamics mimicking the contact of the 
system with a thermal reservoir. These dynamic rules - usually referred to as heat hath 
dynamics - are defined through the transition probabilities 

i j 

where a denotes the actual state of the model and j runs over the nearest neighbors of i. 
In order to verify the coincidence of the stationary distribution Poo{ct) and the equilibrium 
ensemble of the Ising model, it is sufficient to prove that the dynamic processes are ergodic 
and obey detailed balance 

Peqicr)Pa~>a' = Peqicr')Pa'^a ■ (60) 

Detailed balance means that the probability currents between two states are exactly equal 
in both directions, i.e., the currents cancel each other in the stationary state. Heat bath 
dynamics is only one out of infinitely many dynamic processes generating the Ising equilib- 
rium ensemble. Examples include Glauber, Metropolis, and Kawasaki dynamics, as well 
as the Swendsen-Wang and Wolf cluster algorithms. Although these stochastic models 
have very different dynamic properties, they all evolve towards the same stationary state 
which is just the equilibrium state of the Ising model. 
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2.11 Matrix product states 



For the majority of reaction-diffusion models it is quite difficult or even impossible to solve 
the master equation analytically. In some cases, however, it is still possible to compute 
the stationary state of the system. In recent years a powerful algebraic approach has 
been developed by which n-point correlation functions of certain nonequilibrium systems 
can be computed exactly (see Ref. [29| for a general review). This approach generalizes 
states with product measure to so-called matrix product states (MPS) by replacing real- 
valued probabilities with non-commutative operators. Representing these operators in 
terms of matrices it is possible to compute correlation functions by evaluating certain 
matrix products. 



MPS for the exclusion process 

In order to introduce the matrix product technique, let us consider the partially asymmet- 
ric exclusion process in one spatial dimension with asynchronous dynamics and particle 
adsorption (desorption) at the left (right) boundary. The Liouville operator of this system 
is given by 



N~l 



(61) 



i=l 



where Ci describes the hopping in the bulk while Si and Sm are the surface contributions 
for adsorption and desorption, respectively. In the standard basis (see Appendix |A|), these 
operators read 



/O 0\ 

q'^ -q 

g-i q 

0/ 







a 
-a 



-a 



(62) 



The partially asymmetic exclusion process with particle input and output is particularly 
interesting in presence of a current, that is, for q ^ 1- Depending on a and /?, the system 
is in a phase of low, high, or maximal density (see Fig. |5|). It seems to be quite surprising 
that this simple model could be solved only a few years ago by recursion techniques |^,75| 
and, at the same time, by using the matrix product method ||57|,^. 

The solution is particularly simple along the coexistence line between the high and low 
density phases a + P = q — q^^. Here the stationary state \Ps) may be written as 







(63) 



where e = 1/a and d = This state has a product measure, i.e., all spatial corre- 

lations vanish. The constant Z = {e + d)^ normalizes the probability distribution. The 
stationarity C\Ps) = can be verified by proving the relations 



1 



>7V 



1 



d) + id) ^ (-1 



(64) 
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which provide the following cancellation mechanism: First the action of Si generates a 
factor ( Ji) at the leftmost position in the product (|63|). This factor is then commuted to 
the right by successive action of £i, . . . ,Cn-i until it is adsorbed at the right boundary 
by the action of Sm- This cancellation mechanism proves that the homogeneous product 
state ( |63| ) is stationary along the coexistence line a + P = q — q~^. 

For a+f3 ^ q—q~^ the stationary state has no longer the form of a simple product state. 
However, as will be shown below, it can be expressed as a matrix product state. To this end 
we replace the probabilities e and d in Eq. ( |63|) by noncommutative operators E and D. 
These operators act in an auxiliary space which is different from the configurational vector 
space of the lattice model. Introducing boundary vectors {W\ and \V) with (T^IF) / a 
matrix product state may be written as 



1 
Z 



{W\ 



\v) . 



(65) 



Notice that {W\ and \V) are vectors in the auxiliary space while \Ps) denotes the sta- 
tionary probability distribution in configuration space. By selecting the matrix element 
{W\ . . . products of the operators can be mapped to real- valued probabilities^. The 
normalization constant is given by Z = (l^|C^|y), where C = D + E. The cancellation 
mechanism 
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-1 



(66) 



is the same as in Eq. (|6^), leading to the algebra 

DE=D+E 

and the boundary conditions 

{W\E = a-^{W\, D\V) = [3-^\V) . 



(67) 



(68) 



This ansatz is not restricted to the exclusion process but can be applied to any reaction- 
diffusion process with random-sequential updates. It converts the dynamic rules of the 
model into a set of algebraic relations and boundary conditions. The problem of calculating 
the stationary state is then shifted to the problem of finding a matrix representation of 
the algebra. In the present case the quadratic algebra (|6^) can be mapped onto a bosonic 
algebra for which an infinite-dimensional matrix representation exists . For particular 
values of a and j3, however, there are also finite-dimensional representations [^, 78|. For 
example, along the coexistence line a + j3 = q — q~^ the product state (|63|) is nothing 



but a one-dimensional representation of the algebra. Moreover, it can be shown that for 
a + (3 + a(3q = q + q~^ a two-dimensional matrix representation is given by 



E 



1 



\-aq 



D 




{W\ = (1 0) 



1^) 



(69) 



''The matrix product technique can also be applied to models with periodic boundary conditions. In 
this case the matrix elements {W\ . . . \V) have to be replaced by a trace operation Tr[. . . ]. 
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Fig. 8: Exact stationary density profile for the partially asymmetric exclusion process with particle 
adsorption and desorption at the boundaries for q = 2 and A'' — 100 in the low density phase a = 0.45, 
along the coexistence line a = 0.5, and in the high density phase a = 0.55. 



For a given matrix representation the stationary particle concentration profile pf""^ can be 
computed by evaluating the matrix product 



„stat 



{W\C'-^DC^-'\V) 
{W\C^\V) 



In the case of the above 2x2 representation we obtain the exact result 
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aP{l+f5q) 



(70) 



(71) 



(72) 



are the eigenvalues of C. The density profile is shown in Fig. ^ for various values of a, 
visualizing the transition between the low and the high density phase. Similarly, one can 
compute two-point correlation functions 



{W\C'-^DC^-'-^DC^-^\V) 
{W\C^\V) 



(73) 



In this expression the matrix C plays the role of a transfer matrix between the sites i and j. 
Therefore, the long-distance behavior of correlation functions in the bulk will be governed 
by the largest eigenvalue of C. In particular, for any finite-dimensional representation 
of the algebra, the correlation functions in the stationary state will decay exponentially. 
Only infinite-dimensional representations can lead to long-range correlations with power- 
law decay in the stationary state. 

It is also possible to apply the matrix product technique to the totally asymmetric 
exclusion process with sublattice-parallel updates [f7S|-81|, where a different cancellation 
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mechanism is needed ||82|-p5|]. Recently, the matrix product method could even be extended 
to the case of fully parallel updates |4C 



MPS for models with particle reactions 

Most models which have been solved so far by using the matrix product method are 
diffusive systems, i.e., they describe stochastic transport of particles. Up to now only 
one exception is known where particles react with each other, namely the anisotropic 
decoagulation model with closed boundary conditions and random sequential updates. As 
will be explained in the following, this model exhibits a boundary-induced phase transition 
and can be solved by using a generalized matrix product ansatz [^] . 

The anisotropic decoagulation model is defined by the following dynamic rules: 

diffusion: 

coagulation: 

decoagulation: 
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The corresponding Liouville operator reads 








(k + 1) 
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(74) 



q + q 1 / 



Thus, the model is controlled by two parameters, namely the anisotropy q and the decoag- 
ulation rate k. The phase diagram displays two phases, a low-density phase for k < q^ — 1 
and a high-density phase for k > q^ — 1. From the physical point of view these phases are 
different from those observed in the asymmetric exclusion process since here the number 
of particles is not conserved. Notice that the rates for diffusion and coagulation coincide. 
Moreover, all reactions have the same bias. This special choice ensures that the model 
is integrable. Various exact results have been obtained by using IPDF and free-fermion 
techniques ^ . At the critical point Kc = — 1 , the relaxational spectrum 

becomes massless and algebraic long-range correlations can be observed ||5^ 



In order to express the stationary state of the model as a matrix product state of 
the form (|65|), the cancellation mechanism of Eq. ( |66|) has to be generalized. This can 

be achieved by replacing the vector ( j^i ) by ^ ^ ^ j where E and D are two additional 
operators acting in the auxiliary space. The generalized cancellation mechanism reads 
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This ansatz leads to the bulk algebra 
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and the boundary conditions 

{W\E = {W\D = E\V) = D\V) 



0. 



(77) 



A trivial one-dimensional representation of this algebra is given by£' = C = l, E = C = 0, 
describing a system without particles. In the symmetric case q = 1 there is another one- 
dimensional representation E = 1, C = j"^, E = C = 0, corresponding to a homogeneous 



product state with particle density p = k/{1 + k). For g / 1 and /t / 
four-dimensional representation 
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Using this representation, it is easy to compute the stationary particle density 



P 



stat 



(79) 



which is in agreement with the result obtained by using the IPDF method [55|. For 
K 7^ g2 _ X a similar four-dimensional matrix representation can be constructed. 

In the thermodynamic limit the anisotropic coagulation model exhibits a first-order 
phase transition (see Fig. ^). If the decoagulation rate n is small enough, the particles 
are swept towards one of the boundaries where they coagulate. The stationary particle 
density is therefore zero in the thermodynamic limit. Increasing k this region grows until 



its size diverges at a critical value Kc 



1. Above Kc the decoagulation process is 
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strong enough to maintain a non-vanishing density of particles in the bulk. It should be 
emphasized that this type of phase transition is induced by the boundaries. In particular, 
there is no such transition if periodic boundary conditions are used. 

The matrix product technique has also been applied to various other systems such 
as valence-bond-state models ||8^, spin-one quantum antiferromagnets [39, 9C], hard-core 
diffusion of oppositely charged particles systems with fixed [pi] or moving impuri- 



ties [|93|,|9j, as well as n-state diffusion processes [^,^]. Furthermore, a dynamic matrix 
product ansatz has been introduced by which time-dependent properties of the exclusion 
process can be described ^,^|. Although the full range of possible applications is not yet 
known, the matrix product technique seems to be limited to a few classes of models. By 
definition the method is restricted to one-dimensional systems. Moreover, there seems to 
be a subtle connection between matrix product states and integrability. This connection 
is not yet fully understood. In Ref. [|9^ it was shown that the stationary state of any 
reaction-diffusion model can be expressed in terms of a MPS. However, in this generic 
case the corresponding matrix representation depends on the system size and is there- 
fore useless from the practical point of view. A systematic classification scheme for matrix 
product states is not yet known. In this context it is interesting to note that the method of 
dynamic density matrix renormalization allows finite-dimensional matrix representation to 
be detected. As shown in Ref. [100[, the existence of a finite-dimensional MPS is indicated 
by the fact that the density matrix has only a finite number of nonvanishing eigenvalues. 
Thus, by scanning the spectrum over a certain range of the systems parameter space, it 
is possible to search systematically for finite-dimensional matrix product representations. 
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3 Directed percolation 



Spreading processes are encountered in many different situations in nature as diverse as 
epidemics [101|, forest fires [102|, and transport in random media [ 103| , 104|. Spreading 



phenomena are usually characterized by two competing processes. For example, in an 
infectious disease the spreading agent (bacteria) may multiply and infect neighboring 
individuals. On the other hand, infected individuals may recover, decreasing the total 
amount of the spreading agent. Depending on the relative rates for infection and recovery, 
two different situations may emerge. If the infection process dominates, the epidemic 
disease will spread over the entire population, approaching a stationary state in which 
infection and recovery balance one another. However, if recovery dominates, the total 
amount of the spreading agent continues to decrease and eventually vanishes. 

Theoretical interest in models for spreading stems mainly from the emerging phase 
transition between survival and extinction. The simplest model exhibiting such a transi- 
tion is directed percolation (DP)^. In DP sites of a lattice can either be active (infected) 
or inactive (healthy). Depending on a parameter controlling the balance between infection 
and recovery, activity may either spread over the entire system or die out after some time. 
In the latter case the system becomes trapped in a completely inactive state, the so-called 
absorbing state of the model. Since the absorbing state can only be reached but not be 
left, it is impossible to obey detailed balance, i.e., DP is a nonequilibrium process. The 
transition between the active and the inactive phase is continuous and characterized by 
universal critical behavior. 

In many respects, the nonequilibrium critical behavior of DP is similar to that of 
equilibrium models. In particular, it is possible to use the concept of scale invariance and to 
identify various critical exponents. As in equilibrium statistical mechanics, these exponents 
allow phase transitions of different lattice models to be categorized into universality classes. 
As we will see below, DP is a fundamental class of nonequilibrium phase transitions, 
playing a similar role as the Ising universality class in equilibrium statistical mechanics. 
Although DP is very robust and easy to simulate, its critical behavior turns out to be 
highly nontrivial. This is what makes DP so fascinating. 

In this Section we give a general introduction to DP, discussing the most important lat- 
tice models, basic scaling concepts, finite-size properties, as well as mean-field approaches. 
We also summarize various approximation techniques such Monte Carlo simulations, series 
expansions, and numerical diagonalization. Introducing basic field-theoretic methods we 
discuss the critical behavior of DP at surfaces, the early-time behavior, the influence of 
fractal initial conditions, persistence probabilities, and the influence of quenched disorder. 
Finally we review possible experimental realizations of DP and discuss the question why 
it is so difficult to verify the critical exponents in experiments. 



^ For further reading in this field we recommend the review on directed percolation by Kinzel [p5, 105|, 
a summary of ope n problems by Grassberger [ IOC ] , and the relevant chapters in a recent book by Marro 
and Dickman 107 1. 
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isotropic bond percolation directed bond percolation 



Fig. 10; Isotropic and directed bond percolation on a diagonal square lattice with free boundary conditions. 
Open (closed) bonds are represented by solid (dotted) lines. In both cases a cluster, indicated by bold 
hues, is generated from the lattice site in the center. In the directed case the spreading agent is restricted 
to follow sense of the arrows, leading to a directed cluster of connected sites. 



3.1 Directed percolation as a spreading process 
Prom isotropic to directed percolation 

Although DP is often regarded as a dynamic process, it was originally defined as a geo- 
metric model for connectivity in random media which generalizes isotropic (undirected) 



percolation [108, 109 1. Such a random medium could be a porous rock in which neighbor- 
ing pores are connected by channels of varying permeability. An important question in 
geology would be how deep the water can penetrate into the rock. 

In ordinary percolation the water propagates isotropically in all directions of space. 
One of the simplest models for isotropic percolation is bond percolation on a d-dimensional 



square lattice, as shown in the left part of Fig. |10|. In this model the channels of the porous 
medium are represented by bonds between adjacent sites of a square lattice which are open 
with probability p and otherwise closed]^. For simplicity it is assumed that the states of 
different bonds are uncorrelated. Clearly, if p is sufficiently large, the water will percolate 
through the medium over arbitrarily long distances. However, if p is small enough, the 
penetration depth is expected to be finite so that large volumes of the material will be 
impermeable. Both regimes are separated by a continuous phase transition. 



The left part in Fig. |10| shows a typical configuration of open and closed bonds in 
two dimensions. As can be seen, each site generates a certain cluster of connected sites 
corresponding to the maximal spreading range if the water was injected into a single pore. 
The site from where the cluster is generated is called the origin of the cluster. The size 
(or mass) of a cluster is defined as the number of connected sites. Notice that different 
origins may generate the same cluster. Consequently the whole lattice decomposes into a 
set of disjoint clusters. 



Alternatively, we could have blocked sites instead of bonds with a certain probability. The resulting 
model, called site percolation, exhibits the same type of universal critical behavior at the transition. 
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Directed percolation, introduced in 1957 by Broadbent and Hammersley []110| , is an 



anisotropic variant of isotropic percolation. As shown in the right part of Fig. 10, this 
variant introduces a specific direction in space. The channels (bonds) function as 'valves' 
in a way that the spreading agent can only percolate along the given direction, as indicated 
by the arrows. For example, we may think of a porous medium in a gravitational field 
that forces the water to propagate downwards^. Thus, filling in the spreading agent at a 
particular site, the resulting cluster of wet sites is a subset of the corresponding cluster in 
the isotropic case (see Fig. |l^. Since in DP each site generates an individual cluster, a 
decomposition of the lattice into disjoint clusters is no longer possible. As in the case of 
isotropic percolation, DP exhibits a continuous phase transition. 

The phase transitions of isotropic and directed percolation are similar in many re- 
spects. They both can be characterized by an order parameter P^o which is defined as the 
probability that a randomly selected site generates an infinite cluster. If Pqo is finite, the 
spreading agent is able to percolate over arbitrarily long distances wherefore the system 
is said to be in the wet phase. If Poo vanishes, the system is in the so-called dry phase 
where the spreading range is finite. Both isotropic and directed percolation are trivial in 
one dimension: Since an infinite cluster on a line requires all bonds to be open, the wet 
phase consists of a single point p = 1. Another trivial case is the limit of infinitely many 
dimensions. Since each site is connected with infinitely many neighbors, an infinite cluster 
will be generated for any p > 0. Consequently the inactive phase consists of single point 
in phase space, namely p = 0. In finite dimensions 2 < d < oo there is a continuous phase 
transition separating the wet phase from the dry phase at some critical value < pc < 1- 
In the supercritical phase p > Pc the medium is permeable (Poo > 0) while in the subcrit- 
ical phase p < pc the medium becomes impermeable (Poo = 0). As expected, the critical 
threshold pc for directed percolation is larger than in the isotropic case. 

Although isotropic and directed percolation have several common features, their critical 
behavior near the phase transition turns out to be different. In the isotropic case, the 
critical properties are in all directions the same (apart from lattice effects which are usually 
irrelevant on large scales) and hence the emerging long-range correlations are rotationally 
invariant. Because of a duality symmetry, the critical point of isotropic bond percolation 



is Pc = 1/2 |10£]. Moreover, the critical exponents are given by simple rational numbers. 
Contrarily, the critical properties of DP reflect the anisotropy in space, leading to different 
critical exponents. In contrast to the isotropic case, the numerical values of the critical 
point and the exponents of DP are not yet known analytically and seem to be given by 



irrational numbers (see Sec. 3.4). 



Interpretation of directed percolation as a dynamic process 

Regarding the given direction as 'time', directed percolation may be interpreted as a d+1- 
dimensional dynamic process that describes the spreading of some non-conserved agent|^. 
For example, as already mentioned before, DP may be viewed as a simple model for 



epidemic spreading of some infectious disease without immunization [111|. In recent years 
the dynamic interpretation has become increasingly popular, partly because the time- 
dependent formulation is the natural realization of DP on a computer. In what follows 



'^This assumption is highly idealized since water is a conserved quantity. Moreover, the water can even 



flow against the gravity field (see Sec. 3.!:) 



Note that D P d ifi^ers from 'dynamic percolation' which is used as a epidemic processes with immu- 



nization (see Sec. 3.8 
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initial configuration 



final configuration 



Fig. 11: Directed bond percolation in 1+1 dimensions interpreted as a time-dependent stochastic process. 
Open (closed) bonds are indicated by solid (dashed) lines. Filled (hollow) circles denote active (inactive) 
sites. The configuration of the horizontal row at t = is the initial state. Starting from a fully occupied 
initial state the model 'evolves' through intermediate configurations according to the dynamic rules of 
Eq. (KOI) and reaches a final state at t = 3. 



we will adopt the dynamic interpretation. However, one should keep in mind that the 
geometric definition in terms of directed paths is fully equivalent. 



The dynamic interpretation of DP is illustrated in Fig. |11|, where the lattice sites of 
a (l-l-l)-dimensional directed bond percolation process are enumerated horizontally by a 
spatial coordinate i and vertically by a discrete time variable t. Since we use a diagonal 
square lattice, odd and even time steps have to be distinguished. A local binary variable 
Si{t) is attached to each site. Sj = 1 means that the site is active (occupied, wet) while 
Si = denotes an inactive (empty, dry) site. The set s = {si} at a given time t specifies 
the configuration of the system. 

For a given configuration at time t, the next configuration at time t + 1 can be de- 
termined as follows. For each pair of bonds between the sites (i it l,t) and {i,t + 1) 
two random number zf^ G (0, 1) are generated. A bond is considered to be open (with 
probability p) if zf^ < p, leading to the update rule 

{1 if Si-_i(t) = 1 and z^r <c p , 
1 if Si+i{t) = 1 and z+ < p , (80) 
otherwise . 

Thus, directed bond percolation can be considered as a Markov process with parallel 
dynamics. As in any dynamic system, we have to specify the initial state. Common initial 
states are the fully occupied lattice, random initial conditions, and configurations with a 
single particle at the origin (also called 'active seed'). 

Even very simple numerical simulations demonstrate that the temporal evolution of a 
DP process changes significantly at the phase transition. Typical space-time histories for 
random initial conditions are shown in the upper part of Fig. For p < Pc the number 
of particles decreases exponentially until the system reaches the absorbing state, whereas 
for p > Pc the average particle number saturates at some constant value. At the critical 
point the mean particle number decays very slowly and the emerging clusters of active 
sites remind of fractal structures. A similar behavior can be observed if the DP process 
starts from a single seed (see lower part of Fig. |l^). For p < pc the average number of 
particles first grows for a short time and then decays exponentially. For p > pc there is a 
finite probability that the resulting cluster is infinite. In this case activity spreads within 
a certain triangular region, the so-called spreading cone. At p = pc a critical cluster is 
generated from a single seed, whose scaling properties will be discussed in Sec. 



It is often helpful to regard DP as a reaction-diffusion process of interacting particles. 
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Fig. 12: Directed bond percolation in 1+1 dimensions starting from random initial conditions (top) and 
from a single active seed (bottom). Each horizontal row of pixels represents four updates. As can be seen, 
critical DP is a reaction-limited process. 



Associating active sites with particles A and inactive sites with vacancies 0, a DP process 
corresponds to the reaction-diffusion scheme 

self-destruction: A , 

diffusion: + A ^ A + , , , 

offspring production: A 2A , 

coagulation: 2A ^ A . 

To understand this reaction-diffusion scheme, let us again consider the example of directed 
bond percolation. Depending on the configuration of the bonds, each active site (particle) 
may activate two neighboring sites of the subsequent row (next time step) . If both bonds 
are closed, the particle self-destructs. If only one bond is open, the particle will diffuse 
to the left or to the right with equal probability, whereas an offspring is produced when 
both bonds are open. On the other hand, if two particles reach the same target site, they 
coalesce into a single particle, giving rise to the reaction 2A A. This process limits the 
maximal density of active sites. In fact, as will be shown below, the coagulation process is 
the essential nonlinear ingredient of DP. In 'fermionic' models with an exclusion principle 
it is automatically included. However, in 'bosonic' models allowing for an infinite number 
of particles per site one would have to add this process explicitly. 

3.2 Lattice models for directed percolation 

In the literature there is a vast variety of DP models following the spirit of the above 
reaction-diffusion scheme. As we will see below, they all exhibit the same type of critical 
behavior at the transition. The common feature of all these models is the existence of an 
absorbing state, i.e., a configuration that the model can reach but from where it cannot 
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Fig. 13: Transition probabilities in the (l+l)-dimensional Domany-Kinzel cellular automaton. 

escape. In most cases, the absorbing state is just the empty lattice. The existence of 
an absorbing state implies that certain microscopic processes are forbidden (for example, 
spontaneous creation of particles — > ^). In the sequel we will discuss three examples, 
namely the Domany-Kinzel cellular automaton, the contact process, and the Ziff-Gulari- 
Barshad model for heterogeneous catalysis. 

The Domany-Kinzel cellular automaton 

Cellular automata are algorithms that map a configuration of a lattice onto a new configu- 
ration. The automaton evolves in time by iteration of the map. Thus, the time variable t is 
discrete. Usually the map can be decomposed into independent local updates. Since these 
updates can be processed simultaneously, cellular automata can efficiently be implemented 
on parallel computers. Depending on the type of updates, we distinguish between deter- 
ministic and stochastic cellular automata. A general classification of stochastic cellular 



automata was presented by Wolfram [112]. 



Various stochastic cellular automata are known to exhibit a DP transition from a 
fluctuating phase into an absorbing state. One of the simplest models in this class is 
the (l+l)-dimensional Domany-Kinzel (DK) model 105 , 113t| . It is defined on a diagonal 



square lattice and evolves by parallel updates according to certain conditional transition 
probabilities P{si{t + l)|si_i(t), These probabilities depend on two parameters 

and are defined by 

P[1|0,0] = , 

P[1|0,1] =P[1|1,0] =pi , (82) 

where P[0|-,-] = 1 — P[l|-,-]. The corresponding update scheme may be realized by the 
following algorithm (see Fig. |l3|): For each site i we generate a uniformly distributed 
random number Zi{t) S (0, 1) and set 



1 if Si_i(t) / Si+i(t) and Zi{t) < pi , 
1 if Si_i(t) = Si+i(t) = 1 and Zi{t) < p2 , (83) 
otherwise . 



Si{t + 1) = < 

In contrast to directed bond percolation, the DK model depends on two percolation prob- 



abilities pi and p2- The corresponding phase diagram is shown in Fig. 14. It comprises 
an active and an inactive phase, separated by a phase transition line (the solid line in the 
figure). In the active phase a fluctuating steady state exists on the infinite lattice, whereas 
in the inactive phase the model always reaches the absorbing state. The DK model in- 
cludes three special cases. The previously discussed case of directed bond percolation 
corresponds to the choice pi = p and p2 = p{2 — p). Another special case is directed site 
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Fig. 14: Phase diagram of the Domany-Kinzel model. 
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site DP 


0.705489(4) 


0.705489(4) 
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bond DP 


0.6447001(1) 


0.8737620(2) 
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compact DP 


1/2 


1 


M 



Table 1: Special transition points in the (l+l)-dimensional Domany-Kinzel model. 



percolation |^], corresponding to the choice pi = P2 = P- The third special case p2 = is 
equivalent to the rule 'W18' of Wolfram's classification scheme [|112|] . Numerical estimates 
for the corresponding critical points are summarized in Table 0. 

There is strong numerical evidence that the critical behavior along the whole phase 
transition line (except for its upper terminal point) is that of DP. This means that all 
these transition points exhibit the same type of long-range correlations. The short-range 
correlations, however, are non-universal and may change when moving along the phase 
transition line. In order to understand the significance of short-range correlations from 
the physical point of view, let us consider spatial configurations (snapshots) of a critical 
DP cluster. Such configurations are typically characterized by localized spots of activity 
separated by large voids in between. Approaching the phase transition line the average 
size of inactive voids diverges, whereas the mean size of active spots remains finite and 
converges to a certain value {Sact)- By moving along the transition line of the DK model 



the asymptotic average size (Sact) of active spots varies. As shown in Fig. 15, it is minimal 
for p2 = 0, grows monotonically with p2, and finally diverges at the terminal point p2 = 1, 
where the system crosses over to a different type of critical behavior. Thus, by moving 
along the phase transition line the nonuniversal short-range properties change while the 
long-range properties remain unaffected. 

The exceptional behavior at the upper terminal point of the phase transition line is 
due to an additional symmetry between active and inactive sites along the line P2 = 1 (3§1 • 
Here the DK model has two symmetric absorbing states, namely the empty and the fully 
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Fig. 15: Numerical estimates for the average size of active spots {Sact) in the Domany-Kinzel model 
measured along the phase transition line. The insets show typical clusters for three special cases discussed 
in the text. 



occupied lattice. As the corresponding symmetry transformation maps pi to 1 — pi, the 
phase transition line must end in the terminal point {pi,P2) = (1/2,1). As shown in the 



corresponding inset of Fig. |15|, the resulting clusters of active sites are compact. There- 
fore, this special case is referred to as compact directed percolation (CDP) [117| . Unfor- 
tunately this expression is misleading since CDP stands for a universality class which is 
completely different from DP. As can be verified easily, the DK model at the terminal 
point is equivalent to the (l-l-l)-dimensional voter model |Q or the Glauber-Ising model 
at zero temperature. Alternatively, one may describe CDP as a pair-annihilation process 



of diffusing kinks separating inactive and active domains (cf. Sec. 2.4). We may there^ 



fore identify the critical behavior of CDP with the exactly solvable universality class of 



diffusing-annihilating random walks [H. We will come back to CDP in Sec. 3.8 



In more than one spatial dimension the DK model may be defined by local updates with 
conditional probabilities P{si{t + l)|rej(t)) depending on the number nj(i) = Sj{t) 
of active neighboring sites. Thus, the model is controlled by 2d parameters pi, ■ ■ ■ ,P2d- 

P[1|0] = , 

P[l\n]=pn. (l<n<2d) (84) 

Notice that for d = 1 this definition is compatible with the usual definition of the DK 
model in Eq. (^2|). The special case of directed bond percolation corresponds to the choice 
p„ = 1 — (1 — p)" while for equal parameters Pn = P one obtains directed site percolation 
in d+1 dimensions. 

The contact process 

Another important lattice model for DP is the contact process. In contrast to cellular 
automata this model uses asynchronous updates. The contact process was first introduced 
by Harris |111| ] as a model for epidemic spreading without immunization (for a review see 
Ref. | 107[ ). Here the lattice sites represent infected and healthy individuals. Infected 
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Fig. 16: Stochastic processes in the (l+l)-dimensional contact process. Infected sites (bold dots) infect 
their neighbors at rate A/2 and recover at rate 1. The upper part shows the definition of the rules according 
to Eq. (pa). The lower part shows an equivalent definition as a two-site process (see text). 

individuals can either heal themselves or infect their nearest neighbors. Depending on 
the relative rates of infection and recovery, the epidemic disease may either spread over 
the whole population or vanish after some time. In contrast to the DK model, infection 
and healing processes are assumed to occur spontaneously without correlation in space 
and time, i.e. spatially separated processes are not synchronized. To mimic this kind 
of asynchronous dynamics, the contact process uses random sequential instead of parallel 



updates (cf. Sec. 2.2 ) 



The contact process is defined on a d-dimensional square lattice whose sites can be 
either active {si{t) = 1) or inactive (si(t) = 0). For each attempted update a site i 
is selected at random. Depending on its state Si{t) and the number of active neighbors 
ni{t) = X]jG<i> '^i(^) ^ '^^'^ value Si{t+dt) = 0, 1 is assigned according to certain transition 
rates w[si{t) Si{t + dt),ni{t)\. In the standard contact process these rates are defined 
by 

w[0^lM=^n/2d , 

Here the parameter A controls the infection rate and plays the role of the percolation 
probability. For the (l-l-l)-dimensional case the dynamic processes are sketched in the 
upper row of Fig. |l6[ Monte Carlo simulations and series expansions suggest that the phase 



transition in 1+1 dimensions takes place at the critical point Ac — 3.29785(8) [ 118 -12C| 



Whereas computational physicist often prefer the DK model for simulations, the con- 
tact process is more popular in the mathematical community because it is easy to write 
down the corresponding master equation. Following the notation of Sec. |2.2| , the master 
equation of the 1+1 dimensional contact process with periodic boundary conditions is 
given by 



N 



dtPtisi, ... ,SAr) = ^ (2Si - 1) { ASi_i Pt{si, . . . ,Si-2, l,0,Sj+i, ... ,Sn) 

+ ASj+i Pt(si, . . . ,Si_i,0, l,Si+2, . . . ,sn) 
- Pt{si, ... , Si-l, 1, Sj+l, . . . , sn)} , 



i=l 



(86) 



where Pt{si, . . . , sn) denotes the probability to find the system at time t in the configura- 
tion {si, . . . , Sn}. As can be seen, the master equation involves only two-site interactions, 
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as illustrated in the lower row of Fig. |l6[ Using the vector notation of Eq. (P) the corre- 
sponding Liouville operator Lcp = Hi is given by 
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-1 


-1 




1 





1 + A 





-1 


2 








1 + A 


-1 




Vo 


-A 


-A 


2 / 



(87) 



Notice that this operator does not satisfy simple algebraic relations as in Eq. (H), indicat- 
ing that DP is a highly nonintegrable process. Finite-size spectra of Ccp will be analyzed 



in Sec. 3.4. 



The Ziff-Gulari-Barshad model for heterogeneous catalysis 

Many catalytic reactions such as the oxidation of carbon monoxide on a platinum surface 
mimic the reaction scheme of directed percolation. The key property of these reactions 
is the existence of catalytically poisoned states where the system becomes trapped in a 
frozen state. Thus, poisoned states play the role of absorbing configurations. A simple 
model for surface catalysis of the chemical reaction CO + O — > CO2 was introduced in 
1986 by Ziff, Gulari, and Barshad (ZGB) |121U . The model describes a gas composed of 
CO and O2 molecules with fixed concentrations y and 1 — y, respectively, which is brought 
into contact with a catalytic material. The catalytic surface is represented by a square 
lattice whose sites can be either vacant (0), occupied by a CO molecule, or occupied by an 
O atom. The ZGB model evolves by random sequential updates according to the following 
probabilistic rules: 

1. CO molecules fill any vacant site at rate y. 

2. O2 molecules dissociate on the surface into two O atoms and fill pairs of adjacent 
vacant sites at rate 1 — y. 

3. Neighboring CO molecules and O atoms recombine instantaneously to CO2 and 
desorb from the surface, leaving two vacancies behind. 

On the lattice the three processes correspond to the reaction scheme 

CO at rate y , 

+ 0^0 + at rate 1 - y , (88) 
O + CO ^ + at rate 00 . 



Clearly, this reaction is irreversible and thus the dynamic processes do not obey detailed 
balance. Moreover, if the whole lattice is covered either with pure CO or O, the system 
is trapped in a poisoned absorbing state. As shown in Fig. the ZGB model can be 
in three different phases. For y < 7/1 ~ 0.389 the system evolves into the 0-poisoned 
state whereas for y > y2 — 0.525 it always reaches the CO-poisoned state. In between 
the model is catalytically active. The model exhibits two different phase transitions, a 



continuous one at y = yi and a discontinuous one at y = y2- Grinstein et al. [122] 
expected the continuous transition to belong to the DP universality class. In order to 
verify this hypothesis, extensive numerical simulations were performed. Initially it was 
believed that the critical exponents were different from those of DP [123|, while later 
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Fig. 17: Schematic phase diagram of the Ziff-Gulari-Barshad model. The concentrations of oxygen (solid 
line) and carbonmonoxide (dashed line) are plotted versus the CO adsorption rate y. 



the transition at y = yi was found to belong to DP [124|. Very precise estimates of the 
critical exponents were recently obtained in Ref. |125|, confirming the existence of a DP 
transition in the ZGB model. DP exponents were also obtained in a simplified version of 
the ZGB model [ [L26| . However, so far it has been impossible to observe DP exponents in 
experiments. We will come back to this problem in Sec. IsT 



3.3 Phenomenological scaling theory 

In equilibrium statistical physics continuous phase transitions are usually characterized 
by universal scaling laws. For example, the magnetization order parameter in the ordered 
phase of the two-dimensional Ising models vanishes close to the critical point as |T — 
Tc|^, where /3 is a universal exponent. Similarly the correlation length ^, which is the 
characteristic macroscopic length scale of the model, diverges as ^ ~ |T — Tc\~'^. At 
the critical point the correlation length is infinite, i.e., there is no macroscopic length 
scale in the system. As a consequence, the system is invariant under suitable scaling 
transformations. It turns out that a very similar picture emerges in the nonequilibrium 
case. In the following we introduce a phenomenological scaling theory that can be applied 
to DP and other types of phase transitions into absorbing states. 



The critical exponents /?, v^, and i/y 

The order parameter of a spreading process is the density of active sites 

i 

where (...) denotes the ensemble average. Let us first consider the case of an infinite 
system. In the active phase p{t) decays and eventually saturates at some stationary value 
pStat rjij^g stationary density varies continuously with p — Pc and vanishes at the critical 
point (see Fig. |l^) . Close to the transition the order parameter varies according to a power 
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Fig. 18: Stationary density p"*"^' in the active phase of directed bond percolation in 1+1 and 2+1 
dimensions, plotted in a linear and a double-logarithmic representation. 



law 

P^'^'r^ip-Pcf, (90) 

where /? is the critical exponent associated with the particle density. In a double-logarith- 
mic representation the power law behavior manifests itself as a straight line with slope p. 



As can be seen in Fig. |18|, the value of /? depends on the dimensionality of the system. 
The numerical value f3 ~ 0.277 in 1+1 dimensions is comparatively small, indicating a 
significant change of p**"^* near the transition. In 2+1 dimensions a larger value /? ~ 0.58 
is observed. 

In addition, spreading processes are characterized by certain correlation lengths. In 
contrast to equilibrium models without any dynamical aspect, nonequilibrium critical 
phenomena involve 'time' as an additional dimension. Since 'time' and 'space' are different 
in character, we have to distinguish spatial and temporal properties, denoting them by 
the indices + and ||, respectively. In fact, nonequilibrium phase transitions are usually 
characterized by two independent correlation lengths, namely a spatial length scale and 
a temporal length scale ^y. Close to the transition, these length scales are expected to 
diverge as 

e±~|p-Pcr"^ ^h-Ip-pcI-'h (91) 

with generally different critical exponents z^j_ and z^y. In the scaling regime the two corre- 
lation lengths are related by ~ where z = v^jvi, is the so-called dynamic exponent. 
In many models the triplet v^^ is the fundamental set of bulk exponents that 
labels the universality class. Other critical exponents are usually related to these three 
exponents by simple scaling relations (see below). Nonequilibrium phase transitions in 
different physical systems are believed to belong to the same universality class if their 
critical exponents coincide^. In fact, the DK model, the contact process, the ZGB model, 
and a vast variety of other DP models are characterized by the same triplet of exponents. 



Fig. 19 illustrates the physical meaning of the correlation lengths and .^||. As in 



equilibrium statistical mechanics, they are present below and above the critical point. In 
^In addition, it should be proven that universal scaling functions coincide as well. 
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Fig. 19: Interpretation of the correlation lengths and in an almost critical (H-l)-dimensional DP 
process below (left) and above criticality (right). In panels A and B a cluster is grown from a single active 
seed while in panel C a fully occupied lattice is used as initial state. Panel D shows a stationary DP 
process in the active phase. The indicated length scales and must be interpreted as averages over 
many independent realizations. 



the inactive phase clusters originating from a single seed have the typical form of a droplet 
(panel A). Averaging over many independent realizations the lateral size and the lifetime 
of such droplets are proportional to and ^y, respectively. Above criticality the surviving 
clusters grow within a spreading cone (panel B) whose opening angle is determined by the 
ratio The correlation lengths can also be seen if homogeneous initial conditions 

are used. In the inactive phase the scaling length plays the role of a typical decay 
time (panel C), while in the stationary state of the active phase the correlation lengths 
appear as the average sizes of inactive islands (panel D). This interpretation can be easily 
generalized to higher dimensions. 

As suggested by the scaling properties of the density ( |90| ) and the correlation lengths 
(pl|), a spreading process should be invariant under dilatation x Ax accompanied by 
an appropriate rescaling of time and the deviation from criticality A = p — pc- 

x^Ax, t^AH, A^A-i/^'^A, p^A-^/^^p. (92) 

This allows scale-invariant combinations to be constructed such as t/x^, At^/'^H and 
Ax^/^-L. As we will see below, universal scaling functions can only depend on such scale- 
invariant ratios. 



Scaling theory for phase transitions into absorbing states 

So far we have seen that the stationary density in the active phase scales as ~ 
where A = p — Pc denotes the distance from the critical point. A very similar quantity 
is the ultimate survival probability Poo that a randomly chosen site belongs to an infinite 



cluster (cf. Sec. 3J.). In the active phase this probability is finite and scales as 

Poo ~ A^' (93) 

with some critical exponent /?'. Although /? and /3' coincide in the case of DP, they 
may be different in more general contexts, for example, in models with many absorbing 
states. Therefore, phase transitions into absorbing states are generally described by four 
exponents /3, P' , v^, and The different roles of P and (3' become apparent in a field- 
theoretic formulation (see Sec. |3.5| ). It turns out that /? is associated with the particle 
annihilation operator. Therefore, this exponent emerges whenever a particle density is 
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'measured' in some final state. The exponent /?', on the other hand, is associated with 
the creation of particles and thus plays a role whenever particles are 'introduced'. This 
happens, for example, if an initial configuration is specified. In correlation functions, 
which involve creation as well as annihilation operators, both exponents are expected to 
appear. 

Turning to time-dependent scaling properties in an infinitely large system, there are 
two important complementary quantities, namely the particle density p{t) starting from a 
fully occupied lattice, and the survival probability P{t) that a cluster grown from a single 
seed is still active after t time steps. Following the usual scaling concept of equilibrium 
statistical mechanics, both quantities are expected to scale as 

p{t) ^ t-y{A t^l-W ) , Pit) ^ t-'g{A t^l-W ) , (94) 



where a and 5 are certain critical exponents for decay and survival, respectively [ 127| , |128|| . 



/ and g are universal scaling functions, i.e., they have the same functional form in all DP 
models. For small arguments they both tend to a constant, whereas for large arguments 
they scale in a way that the time dependence drops out: 

/(C) 5(C) -C''^". (C-oo) (95) 

In the active phase the two quantities therefore saturate at p^^"'^ = p{oo) ~ A^'^H and 
Poo = ^'(oo) ~ A'^'^il . Comparison with Eqs. (@) and (H) yields 

a = l5/v\\ , 5 = lv\\ . (96) 

An important quantity that combines both creation and annihilation of particles is the 'pair 
connectedness function c(x', i', x, i), which is defined as the probability that the sites (x', t!) 
and (x, t) are connected by a directed path of open bonds. Since the pair connectedness 
function is translationally invariant in space and time, we may also write c(x',t',x, t) = 
c(x — x', f — t'). Starting from an initial condition with a single active site at the origin x' = 
0, the pair connectedness function c(x, t) is just the density of active sites in the resulting 
clusters averaged over many realizations of randomness. Because of scaling invariance, the 



pair connectedness function c(x, t) obeys the scaling form [ 127 | 



c(x, t) ~ t^-'^l^ F{x/t^l' , A t^/'^W ) , (97) 

where d denotes the spatial dimension and z = The so-called critical initial slip 

exponent 6 describes the growth of the average number of particles as a function of time 



(see Sec. |3.6D . In order to determine 6 we note that in the active phase A > surviving 
clusters will create an average density p'^*"* ~ A^ in the interior of the spreading cone (cf. 
panel B of Fig. [l9| ). Thus the autocorrelation function c(0,t) should saturate at the value 

c(0,oo) = lim c(0,t) ~ A^+z^' . (98) 

t— ►oo 

On the other hand, the scaling form ( [97|) implies that c(0,t) saturates in the active phase 
at a constant with the scaling behavioifn 



c 



(0,oo) ~ A^ll('^/^-''). (99) 



"'To prove this relation, notice that F(0,C) ~ ^"(^"''/'^'''ll for large values of C- 
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Comparing the two expressions we obtain the generalized hyperscaling relation 129 | for 
phase transitions into absorbing states 

= (100) 

z z^ii 

It should be noted that the scaUng argument (^) rehes on the assumption that the cluster 
spreads around the origin, i.e., the spreading cone surrounds the origin. In sufficiently high 
spatial dimensions, however, the cone becomes sparse and diffuses away from the origin 
so that the autocorrelation function c(0, oo) vanishes. For example, in a DP process this 
happens above the upper critical dimension dc = 4. In fact, the generalized hyperscaling 
relation ( |100 ) turns out to be valid only below the upper critical dimension of the spreading 



process under consideration. 

The scaling theory outlined above assumes the system size to be infinite. For finite 
systems sizes the scaling functions also depend on the invariant ratio (,'j_/N = f^^^/N, 
where N = L"^ is the total number of sites. The generalized scaling forms read 

p{t) ~ t-^/'^it /(At^/^ll, t^/VAf) , (101) 
P{t) ~ t~^>ii5(Ati/^ll, t'^/ViV) , (102) 
c(x,i) ~ t-^^+f^'^/^w F{x/t^/',At^/''\\,t'^/yN) . (103) 

Derived scaling properties 

The scaling behavior of various other quantities can be derived directly from the scaling 
relations (Il0l| )- (|103D . For example, the mean cluster mass M is given by the total integral 
of the pair connectedness function 



M = I d'^x I dtc{x,t) . (104) 



oo 







Inserting the scaling relation (|97| ) and substituting the scaling variables we obtain a scaling 
law for the average cluster mass measured in an infinite system below criticality: 

M / d'^x / d^^^-^/^F(x/^^/^A^l/'^||) ~ |Ar''ii(i+^) . (105) 







Similarly, the mean survival time T, the mean spatial volume V, and the mean size S" of a 
cluster in the inactive phase are given by 

T = J dtP{t) = J dtt-^ GiAt^/^w) ~ lAI-^iK^--^) , (106) 
V = J dtP{tyi^-^ = j dtf^/'-^-^ GiAt^l^w) ~ ^ (107) 

S = j dtP{tyl^ = j dtf^'^-^ GiAt^/^W) ~ \l^\-n{d/z+l-5) _ (-^Qg^ 

For these quantities, we obtain the following scaling relations: 

(109) 
(110) 
(111) 

(5' . (112) 



M ~ 


|A|-T , 


7 = 


^1 


{i + e) = 


z^ll + dv± 
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y\\ - P' , 
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Spreading processes in an external field 

Let us finally consider a spreading process in an external field h. Using the particle 
interpretation, such a field may be realized by spontaneous creation of particles — > ^ 
at rate h during the temporal evolution. Clearly, spontaneous particle creation destroys 
the absorbing state and therefore the transition itself. That is, the external field drives 
the system away from criticality. For small h the resulting distance from criticality obeys 
certain scaling laws. 

In principle the presence of an external field requires to introduce another independent 
critical exponent for the coupling constant. In the case of DP, however, this exponent is not 
independent, it is rather identical with the mean cluster size exponent 7. To understand 
this relation, let us consider the stationary state of a subcritical DP process in presence of 
a weak field. Obviously, a site can only become active if it is connected with at least one 
other site backwards in time where a particle was spontaneously created by the external 
field. Since the number of such sites is equal to the cluster size, the probability to become 
active is given by p^*""* ~ 1 — (1 — h)^^'^^\ For weak fields, the stationary density is 
therefore linear in h: 

~ /iM(A). (A<0) (113) 
Consequently, the susceptibility of a supercritical DP process scales as 

X = ^P'"''^\^\-^ ■ (114) 
Invariance under rescaling ( p2|) requires the external field to change as 

/i^ A^>^-^-'^/i = A-'^/'^^/i. (115) 
Thus, at criticality the stationary response of a DP process is given by 



More generally, we may extend the scaling forms ( 10lD -( p!03[ ) by including the scale- 



invariant argument /it'^^'^li . For example, the density p(t) evolves as 

pit) ~ r^/'^ii /(A ti/'^ii , f^/'/N, ht^l^W ) . (117) 



Time reversal symmetry of directed percolation 

As shown in Ref. [127|, a special symmetry of DP under time reversal implies the additional 



scaling relation 

/? = /?'. (118) 

This is the reason why DP is characterized by only three instead of four critical exponents. 
In order to understand this duality symmetry from the physical point of view, let us 
consider the special case of directed bond percolation. By reversing the arrows shown 



in the right part of Fig. 10 one obtains a directed bond percolation process that evolves 
'backwards' in time. Obviously, the reversed process follows exactly the same probabilistic 
rules as the original one. Moreover, if two sites were connected by a directed path in the 
original process, they will also be connected in the reversed process. Hence, if the reversed 
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process was started from a fully occupied lattice at t > 0, the resulting active sites at t = 
would be precisely those sites which ~ in the original process - would generate clusters 
that are still alive at time t. Therefore, in the case of directed bond percolation we obtain 

P{t) = p{t), (119) 

i.e., the survival probability of a single seed P{t) is exactly equal to the density of ac- 
tive sites p{t) in a DP process starting with a fully occupied lattice. Thus, in the active 
phase the two quantities saturate at the same value Poo = /o**'^* wherefore the correspond- 
ing critical exponents /3 and /?' have to be identical. It should be emphasized that this 
time reversal symmetry of DP is nontrivial and does not hold for other systems such as 
models with several absorbing states |129(| or spreading processes with a fluctuating back- 
ground [|130| . Together with Eq. ( |118| ) the generalized hyperscaling relation (|100| ) reduces 
to the DP hyperscaling relation 

e = d/z-26 = (120) 

Consequently, the autocorrelation function c(0, t) for DP saturates in the active phase at 
the value c(0,oo) - A^/^. 



The DP conjecture 

One of the most fascinating properties of DP models is their robustness with respect to the 
microscopic dynamic rules. In fact, the DP class covers a wide range of models. It includes, 
for example, the vast majority of spreading models such as the contact process [|, 131], 
epidemic spreading without immunization |132| , and forest fire models |102, 133| , 134]. 



Moreover, the DP class includes models for catalytic reactions J 121, 135, 136 
ing particles ]]137|, as well as branching-annihilating random walks with odd number of 
offspring ]13S- |140|] . Furthermore, certain growth processes 1141, 142] and coupled map 
lattices with asynchronous updates 1143] display DP behavior. In fact, this list is far from 
being complete. 

The variety and robustness of DP models led Janssen and Grassberger to the conjec- 
ture that a model should belong to the DP universality class if the following conditions 
hold p44|,[l45[]: 



1. The model displays a continuous phase transition from a fluctuating active phase 
into a unique absorbing state. 

2. The transition is characterized by a positive one- component order parameter. 

3. The dynamic rules involve only short-range processes. 

4. The system has no special attributes such as additional symmetries or quenched 
randomness. 



Although this conjecture has not yet been proven rigorously, it is highly supported by 
numerical evidence. In fact, DP seems to be even more general and may be identified in 
systems that violate some of the four conditions, for example in certain models with non- 
unique 1|146| " |149| or fluctuating passive states 1 129 |. Even complicated spreading processes 
with several spreading agents and multicomponent order parameters were shown to exhibit 
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DP behavior iT2^ , [T24i|l4^ , [T50| -p^. Some models with infinitely many absorbing states, 
that were initially thought to belong to different universality classes, were later found to be 
in the DP class as well (see Sec. 3.8). Not only the bulk exponents /?, u\\ are universal 
but also other quantities as, for example, scaling functions and higher moments of the 
order parameter |12C]. 

It is remarkable that DP is one of very few critical phenomena in 1+1 dimensions 
which has not yet been solved exactly. Despite of its simplicity and robustness it seems to 
be impossible to compute the critical exponents exactly. In fact, the numerical estimates 
suggest that the critical exponents are given by irrational numbers rather than simple 
rational values. The lack of analytical results may be related to the fact that DP - in 
contrast to ordinary (isotropic) percolation - is not conformally invariant since there is 
no symmetry between 'space' and 'time'. Attempts to replace conformal invariance by an 
anisotropic scaling theory have not yet been successfully applied to DP [154|. 

Only few exceptions of DP are known so far. They all violate at least one of the 
four conditions listed above. For example, a different universality class emerges when 
the system has two or more symmetric absorbing states (see Sec. 4.2). Another example 
is the activated random walk model with a conserved order parameter (see Sec. |4.3D 
Different universal properties are also encountered in models where activity spreads over 
long distances by Levy flights (see Sec. |4.1|). 



3.4 Estimation of the critical exponents 

Directed percolation is one of very few critical phenomena whose critical exponents in 1+1 
dimensions are not known exactly. However, thanks to extensive numerical simulations, 
transfer matrix techniques, series expansions, and field-theoretical calculations the critical 
exponents have been estimated in various dimensions to an extremely high accuracy. This 
subsection briefly summarizes the available methods and the most precise estimates. 

Mean-field approximation 

In order to estimate the critical exponents /3, and v\\ of directed percolation, let us 
first consider a simple mean-field (MF) approximation. Denoting by p{t) the density of 
active sites at time t averaged over the entire system, the MF rate equation for the contact 
process ( |85| ) reads 

dtp{t) = {\-i)p{t)-\p\t). (121) 

This equation has two stationary solutions, namely p''*"* = and p'^^"'^ = (A — 1)/A. Hence 
the mean-field critical point is Ac = 1. The solution p**"* = represents the absorbing 
state from where the system cannot escape. In the inactive phase A < Ac the absorbing 
state is stable while the other solution with negative density is unphysical. In the active 
phase A > Ac the absorbing state becomes unstable against small perturbations while the 
second solution represents a stable active state. Near criticality the stationary density 
vanishes linearly as p**"* ^ \ — Therefore, the mean field density exponent is given 
by (5^^ = 1. On the other hand, Eq. ( |121| ) implies the density to decay in the inactive 
phase asymptotically as p{t) ~ exp(— |A — \c\t) ~ exp(— t/^y), hence ~ |A — Ac|~"^ and 
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In order to determine the spatial scaling exponent v^^, the mean- field rate equa- 
tion ( |12lD has to be extended by a term for particle diffusion 

dtp{x, t) = Ap(x, t) - Ap2(x, t) + ^^(x, t) , (122) 

where D is the diffusion constant and A = A — Ac is the deviation from criticality. In a 
lattice model this term corresponds to nearest-neighbor interactions. According to Eq. 
the density p(x, t) changes under rescaling ( |9^ ) as 



p(x, t) ^ A-^/''^p(Ax, AH) . (123) 
Simple dimensional analysis shows that Eq. (|122| ) is invariant under rescaling if 



The above mean field approximation becomes exact in the limit of infinitely many di- 
mensions. In a finite-dimensional contact process it is not clear whether the mean field 
approximation still applies, it is even not obvious that a continuous phase transition still 
exists. However, Liggett was able to rigorously prove the existence of a phase transi- 
tion for a contact process in d > 1 dimensions. As will be shown below, the mean field 
exponents turn out to be exact for d > 4, where dc = 4 is the upper critical dimension of 
DP. Note that the mean field exponents satisfy the hyperscaling relation ( |120D precisely 
in d = 4 dimensions. 

In order to go beyond the standard mean field approximation in low-dimensional sys- 
tems spatial correlations have to be taken into account. An improved mean field ap- 
proximation for the contact process in 1-1-1 dimensions was developed by Ben-Naim and 
Krapivsky |155| , who expressed the temporal evolution of empty intervals on the lattice by 
an infinite hierarchy of differential equations. This approach is very similar to the IPDF 



technique introduced in Sec. 2.8 Approximating the probability to find pairs of neighbor- 
ing empty intervals by the product of single-interval probabilities, they derived a set of 
equations which can be solved exactly. In this approximation the critical exponents are 
given hy P = iy± = 1, and i/y = z = |. Odor could improve these estimates by using a 



generalized mean field approximation combined with coherent anomaly techniques |156], 
reaching an accuracy of almost 1%. 

Monte Carlo simulations with homogeneous initial conditions 

In order to study nonequilibrium phase transitions quantitatively, numerical techniques 
such as Monte Carlo simulations have become an important tool. Because of the steadily 
growing computer capacity the critical exponents can nowadays be estimated within a few 
per cent, in some cases even up to four digits. Further progress is expected as reaction- 
diffusion models can easily be simulated on parallel computers with a large number of 



simple processors ||45 |. 



The simplest numerical method that allows the critical exponents to be estimated is 
a Monte Carlo (MC) simulation starting with a fully occupied lattice. This technique is 
based on the scaling properties of Eq. (101). In a first series of simulations the critical 



percolation threshold has to be determined by measuring deviations from the asymptotic 
power-law decay p{t) ~ t^^ in a sufficiently large system. To this end p{t) is plotted 
versus t in a double- logarithmic graph (see Fig. |20|a). Positive (negative) curvature for 
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Fig. 20: Ordinary Monte Carlo simulations of a (l+l)-dimensional directed bond percolation process 
starting from a fully occupied lattice. Part (a) shows the particle density as a function of time and 
demonstrates the determination of the critical point. Parts (b) and (c) show data collapses for off-critical 
and finite-size simulations, respectively (see text). 



large t indicates that the system is still in the active (inactive) phase. It should be carefully 
analyzed to what extent the estimate depends on the system size used in the simulation. If 
finite-size effects play a role, extrapolation techniques should be used in order to improve 



the estimate [157|. 



Having determined pc and 5 the exponent P may be estimated by measuring the 
stationary density of active sites p^^""^ ~ A'^ in the active phase. However, this type of 
estimate is known to be quite inaccurate since the equilibration time to reach the stationary 
state grows rapidly as the critical point is approached. This critical slowing down can be 
controlled by plotting p{t) versus t A'^H for different values of A and tuning i/y in a way 
that all curves collapse (see Fig. pO|b). The exponent /? is then given by /3 = Suii. Finally, 



the exponent i'± can be determined by finite-size simulations. According to Eq. (101) 



p{t) has to be plotted against t/N^^'^ for various system sizes (see Fig. 20^). By tuning z 
the data points collapse onto a single curve which gives an estimate for i'± = I'w/z. 

Monte Carlo simulations with localized initial conditions 

More accurate estimates for the critical exponents can be obtained by dynamic simulations 
starting from a single particle (active seed) [|127| . This technique exploits the scaling 
properties of the pair-connectedness function c(x, t). Starting from a single particle, one 
measures the survival probability P{t), the number of active sites N{t), and the mean 
square of spreading from the origin B?{t) averaged over surviving runs. According to 



Eqs. (102) and (|103| ) these quantities obey the scaling forms 

P{t) ~ t-^ g{A t^/^W , t'^l'/N) , (125) 
N{t) = [ d'^xc{x,t) ~ t^F(At^/^li, f^/'/N) , (126) 



R^t) = {x\t)) = J d'^xx2c(x,t) ~ t^/' F{At^/''\\, f^'^/N) . (127) 

At criticality, they are expected to display asymptotic power laws 

p(^)~^-^ iv(^)~^^ R^{t)^fi\ (128) 
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Fig. 21: Dynamic Monte Carlo simulations of a directed bond percolation process. The effective exponents 
S(t), 9{t), and 2/z{t) are shown as a function of time for p — 0.6446, 0.6447, 0.6448. The arrows mark the 
actual values of the critical exponents. 



i.e., they show straight hnes in double logarithmic plots. Off criticality, the lines are 
curved, allowing a precise determination of the percolation threshold pc- Technically it is 
often useful to consider local slopes of these curves by introducing effective exponents 

log^J P(t)/P(t/b)) , ^ 

-6(t) = V w Jj ^29 

logio b 

and similarly 9{t) and 2/z{t), where log^g ^ is the distance used for estimating the slope. 
Plotting the local slopes as functions of 1/t, the curves may be extrapolated to t — > oo, as 
illustrated in Fig. 21. The same method works also in higher dimensional systems [158]. 



In order to improve the estimates, it is useful to eliminate the curvature of the data points 
at criticality by plotting the quantities (128) against 1/t^ instead of 1/t. 



Since the spatial size of the growing cluster at a given time is finite, the simulation 
can be accelerated considerably by storing the coordinates of active particles in a dynam- 
ically generated list. Especially at criticality, where the density of particles is low, such 
algorithms are much more efficient. Moreover, finite-size effects are eliminated completely. 

Numerical diagonalization 

The critical exponents may also be approximated by numerical diagonalization of the 
evolution operator. Although the resulting estimates are usually less accurate than those 
obtained by other methods, this technique is of conceptual interest. Let us, for example, 
consider the (l-l-l)-dimensional contact process on a finite lattice with sites and periodic 
boundary conditions which is defined by the Liouville operator (|87|). Solving the eigenvalue 
problem 

CcpW = \^^) (130) 

we obtain a spectrum of eigenvalues {/^j}, as shown in the left panel of Fig. As in all 
reaction-diffusion models, the lowest eigenvalue //q vanishes. The corresponding stationary 
state I'i/'o) is the absorbing state of the contact process. The other eigenvectors represent 



the relaxational modes of the system. As can be seen in Fig. all of them have a 
short lifetime except for the first excited state \tpi) whose eigenvalue tends to zero as 
A increases. This eigenvector represents the active state of the system. In finite systems 
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Fig. 22: Numerical diagonalization of the Liouville operator of the contact process. Left panel: Low-lying 
part of the relaxational spectrum of the (l+l)-dimensional contact process on a lattice with 8 sites and 
periodic boundary conditions. Right panel: Data collapse according to the scaling form (131). 



there is always a finite probability to reach the absorbing state, hence /xi > 0. In infinite 
systems, however, this eigenvalue decreases with A and vanishes at the critical point. Since 
the amplitude of iV'i) decays in the inactive phase as e~^^^ , we may identify fi^"^ with the 
temporal scaling length ^y. In an infinite system we therefore expect /ii to decrease as 
Hi ~ |A — Acl^'ll for A < Ac and to vanish for A > Ac- The corresponding scaling form reads 

IJLi^ N-'l'^h{/\N^I'^''^) , (131) 

where A = A — Ac- Thus, by plotting fiiN^^'^ against AN^^'^'^-'-, the exponents z and u±_ 
can be determined by data collapse, as demonstrated in the right panel of Fig. In 
order to determine the exponent /3, it would be necessary to analyze the components of 
the eigenvector with respect to the particle density in the active phase. A similar 
analysis of DP models with parallel updates, which are defined by transfer matrices instead 
of Liouville operators, can be found in Ref. [ p.59| . 



Density matrix renormalization group methods 

The method of numerical diagonalization can be improved considerably by using density 
matrix renormalization group (DMRG) techniques. The concept of DMRG was introduced 
in 1992 by White [160| in the context of equilibrium statistical physics as a tool for the 
diagonalization of quantum spin chains. The main idea is to prolongate a given spin chain 
by inserting additional spins and to reduce the resulting configuration space by a suitable 
projection mechanism, keeping only the most relevant eigenstates. This renormalization 
procedure is then repeated many times and the spectrum of the iterated Hamiltonian is 
analyzed. Recently DMRG techniques have also been applied to various (l-l-l)-dimensional 



nonequilibrium systems |161~ |165|| (see |100f| for a general overview). The method yields 
surprisingly accurate results. For example, Carlon et al. |162| ] were able to estimate the 
critical exponents of DP by = 0.249(3), = 1.08(2), and z = 1.580(1), deviating 

from the currently accepted values by less than 1.5%. 



Series expansions 
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The most precise estimates of the DP exponents in 1+1 dimensions have been obtained 



by series expansions |119|. This technique is very similar to low- or high-temperature 
expansions in equilibrium statistical physics. As an example let us consider the (1+1)- 
dimensional contact process. Its Liouville operator (|87|) may be separated into two parts 
£(A) = Cq + A£i, where Cq and Ci describe spontaneous self-destruction A ^ and 
offspring production A 2A, respectively. The basic idea is to regard A as a small 
perturbation and to express physical quantities as power series in A. To this end it is useful 
to introduce the Laplace transform of the probability distribution \Pt) and to expand it 
in powers of A: 

/•oo 

\P{s)) = / dte-^' \Pt) = V A" \Pnis)) . (132) 

•^0 n=0 

By applying C{X) from the left one can easily derive the recursion relation 

where |Po) denotes the initial particle configuration. Hence, if the process started from a 
configuration with a single particle, the vector |P„(s)) describes an ensemble of configura- 
tions with at most n particles. It is therefore possible to explicitly construct the vectors 



\Pn{s)), as described in detail in Ref. [119|. 



The above expansion allows the temporal integral of any observable X{t) = {l\X\Pt) 
to be expressed as a power series in A (for notations see Appendix ^): 

/■oo ^ 

/ dtX{t) = lim {l\X\P{s)) = y A" lim {l\X\Pn{s)) . (134) 

•^^ n=0 

Let us, for example, consider the survival probability P{t) that the system has not yet 
reached the absorbing state at time t [cf. Eq. (94)]. Using the vector formalism this 
quantity may be written as 

Pit) = 1 - {0\Pt) = {l\Pt) - {0\Pt) , (135) 

where (0| = (1,0,0,... ,0) denotes the absorbing state. The critical exponents can be 
estimated as follows. On the one hand, the mean survival time T of clusters in the 
inactive phase can be expanded in powers of A by 

/oo °° / \ 

dtPit) = ^ A" Imif (l|P„(s)) - (0|P„(s)) . (136) 
n=0 '^^ V / 



On the other hand, according to Eq. (110) we have T ~ (-A)^''*"" so that 



-J^ In T ~ ^ ^ + const . (137) 

uA Ac — A 

Therefore, in order to estimate Ac and /? — z^y, three steps have to be taken. At first, the 
vectors |P„(s)) have to be determined by iterating Eq. ( [133D up to order Umax- Although 
this recursion relation is quite complicated, it is still simple enough to be implemented on a 



computer (for example, in Ref. [Ill9|| the iteration was carried out up to order Umax = 24). 
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1 - e/4 - 0.01283 

e/12 + 0.03751 
1 + e/6 + 0.06683 e^ 
1- e/12 + 0.03317 e^ 
2 + eVl8 



Table 2: Estimates for the critical exponents of directed percolation obtained by mean field (MF), 
improved mean field (IMF), numerical, as well as field-theoretical methods. 



Next, one has to express T as a power series in A. Finally, Ac and [5 — v\\ can be estimated 



by determining the location and the amplitude of the singularity in Eq. (137) and by using 
a Pade approximation |^. Since the singularity is approached from the inactive phase, we 
are dealing with a subcritical expansion. Similarly one may also consider the supercritical 
case by expanding = /u£o + -^i in powers of fi. 

A general review on series expansion can be found in Ref. 0]. Series expansions were 



applied to (l-l-l)-dimensional DP firstly in Ref. 164 |, where the critical exponents could 
be determined with a relative accuracy of about 10~^. Refined simulations Jl65| ] led the 
authors to the conjecture that the DP exponents should be given by the rational values 
P = 199/720, z^_L = 26/15, and u\\ = 79/72. In a sequence of papers (ll|,|ll|,[l66|,|l6| the 
error margins could be further reduced down to 10~^ . . . 10~^. These improved estimates 
showed that the conjectured rational values were incorrect, indicating that the critical 
exponents of DP could be given by irrational numbers. This should be taken as a warning 
that critical exponents of non-integrable systems are usually not given by simple rational 



values. Currently, the most precise estimates are given in Ref. |168]. Series expansions for 



DP were also performed in two spatial dimensions [169|. In addition, the exponents were 
found to be independent of the type of lattice under consideration. For easy reference we 
listed the most precise estimates in Table |2|. 



Field-theoretical approximations 

By a field-theoretical renormalization group calculation (see Sec. |3.5| ) it is possible to com- 
pute fiuctuation corrections of the critical exponents close to the upper critical dimension 

dc 



4 in powers of e 
are given by 



dc — d. In a two-loop approximation |144, 171 1 these corrections 



53, 4 
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2+^/16 + 
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6/12)2 + 0(63), 
(6/12)2 +0(6^), 
(6/12)2 +0(6^). 



(138) 



For d <2 these approximations are quite inaccurate. However, in three spatial dimensions, 
where numerical simulations and series expansions are difficult to perform, the two-loop 
approximations are regarded as the most precise estimates available. 
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3.5 Field-theoretic formulation of directed percolation 



The robustness of the DP universahty class can be partly understood by studying the 
corresponding field theory. It is interesting to note that the field theory of DP was first 
discovered in a quite different field of physics, namely in the context of hadronic interac- 
tions at ultra-relativistic energies. In order to predict the cross sections of such particles 
at high energies quantitatively, a field-theoretic approach, called Reggeon field theory, was 
developed in the 70's (see Refs. jl72~177|, a general review is given in [I78|). Surprisingly 
it took almost another ten years to realize that Reggeon field theory was nothing but a 
field-theoretic realization of the contact process 1 179 -181], sometimes also called Gribov 
process | 182 , 183]. In the following we sketch the main ideas of a field-theoretic approach 
to DP. 



The DP Langevin equation 

The Langevin equation of motion for directed percolation can be derived directly from the 



master equation for the contact process [144] and reads 

dtpix,t) = Kp(x,t)-A/92(x,t) + Z)VV(x,t) + C(x,t) . 



(139) 



It differs from the mean field equation ( 122 ) by a density-dependent Gaussian noise field 
C(x, t), which is defined by its correlations 



(C(x,t)) = 0, 
(C(x, t)C(x', t')) = T p(x, t) d^i^ - x') 6{t - t') . 



(140) 



Since the amplitude of C(x, t) is proportional to ^Jf){^^, the noise is said to be multiplica- 
tive. This ensures that the absorbing state /?(x, t) = does not fluctuate. The square-root 
behavior stems from the definition of p(x, t) as a coarse-grained density of active sites 
averaged over some mesoscopic box size. Only active sites in this box give rise to fluctua- 
tions of the density, generating a bounded uncorrelated noise. The noise field C(x, i) can 
be viewed as the sum of all these noise contributions in the box. According to the central 
limit theorem, if the number of particles in the box is sufficiently high, C(x, t) tends to a 
Gaussian distribution with an amplitude proportional to the square root of the number 
of active sites in the box. This type of noise has to be distinguished from other nonequi- 
librium systems with multiplicative noise where the noise amplitude is proportional to 



the field /o(x, i) itself without square root ]184, 185]. These systems do not belong to the 
DP class, rather they are related to the KPZ universality class [186]. The DP Langevin 



equation was also tested numerically in Ref. ]187], confirming that the critical exponents 
are in agreement with those of ordinary DP lattice models. Note that in contrast to the 
annihilation process discussed in Sec. |2^, the noise (|140| ) is real due to positive density 
correlations in the bulk. 



The Langevin equation ( |139D can be seen as a minimal equation needed to describe 
DP. It may also include higher order terms such as /9^(x,t) or V^/9(x, i), but these contri- 
butions turn out to be irrelevant under renormalization group transformations. The same 
applies to higher-order contributions to the noise. These additional terms account for 
(nonuniversal) short-range correlations while they are irrelevant on large scales. In fact, 
the robustness of DP originates in the irrelevance of higher-order terms in the Langevin 
equation. 
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Relation to Reggeon field theory 

In field-theoretic calculations it is often more convenient to characterize the dynamic 
system by a partition sum Z. The sum is carried out over all realizations of the field 
(j){x,t) = p(x, and the noise C{x,t), weighted by an appropriate effective action. More 
precisely, the partition sum is defined as the integral over all realizations of the field (j){x,t) 
and the noise t) which satisfy the Langevin equation. Therefore, we may write the 
integrand as a 5- function with Eq. (|139|) as its argument: 



Z ~ j DCP[C] j D(f) /[(/>] 6(^dt(l)- DV^(P - + -(^y (141) 

Here D( and Dcf) denote functional integration, P[C] is the probability distribution of the 
noise field, and /[(/>] stands for an appropriate Jacobian which turns out to be irrelevant 
in the present problem. As shown in Appendix ^ it is possible to integrate the noise by 
introducing a Martin-Siggia- Rosen response field 0(x, t). The resulting action S = So+Sint 
with 



d^x dt 0(x, t) rdt - DV"^ - k 0(x, t) 



^ I (fx dt </>(x, t) [ 0(x, t) - </>(x, t) ] (/.(X, t) 



(142) 
(143) 



is the effective action of Reggeon field theory [178|. In momentum space it may also be 
written as 



5*0 

Sint [4>, ' 



dki^ 



-k, —uj) {—iTuj + Dk^ — K)(f){k, Lo) 



dkoj / dt 



-k, —Lo) 4>{—k\ —uj') 



+ k' , io + uj') — <p{k + k' , laJ + io') 



(144) 
(145) 



where dk^j = (2'7r) d kdto. Formally the free part of the action can be expressed as 

~. 1 



Sn 



So{k,uj) 



dkw '^{-k,-uj)So{k,uj)^{k,oj) , 

Dk'^ - K + iruj 

Dk"^ — K — iroj 



(146) 



where $ = (</>,(/)). Introducing external currents J{k,uj) and J{k,uj) we can rewrite the 
partition sum as 



Z[J,J]~ D<PD4>I'[^J] exp(-So[(P,4>]-S^nt['l>,4>]+ d'^x dt iJ<p + J4^) 



exp( —C 



-int 



-_5_ _5_ 



) ^^^(2 / '^kwJ'{-k,-uj)QQ{k,uj)J{k,uj)^ , (147) 



where the free propagator Qq = Sq is given by 



Go{k,uj) 



Go{k,uj) 
Go{-k,-Lo) 



(148) 
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Fig. 23: Critical DP cluster (grey) and the backbone of the pair connectedness function c(x', t', x, t), 
illustrating the physical meaning of Feynman diagrams in directed percolation. 

with Go{k,u}) = {Dk"^ — k — iTLo)~^ . Because oiQo{q,uj) ^ Q^i^q^—iS) DP is an irreversible 
process. 

Cluster backbone and Feynman diagrams 

Before turning to field-theoretic renormalization group techniques let discuss the physical 
meaning of Feynman diagrams in directed percolation. The full propagator of the field 
theory is the pair connectedness function c(x',i',x, i) which is defined as the probability 
that two sites (x',t') and (x, are connected by a directed path of open bonds (see 
Sec. p. 31 ). In a given realization of open and closed bonds there may be several possible 
directed paths connecting the two sites, as illustrated in Fig. |2^. The union of all possible 
paths constitutes the so-called backbone of the pair connectedness function |188| . More 
precisely, the backbone consist of all sites that are connected with the sites (x', t') and (x, t) 
by a directed walk, i.e., we cut off all dangling ends of the cluster. From the topological 
point of view the backbone is a directed graph consisting of branching and merging lines. 
Because of the duality symmetry ( |lig| ), it is statistically invariant under time reversal. 

In principle the full propagator is given by a weighted sum over all possible backbone 
configurations. Fluctuation effects are mainly due to the influence of closed loops. Above 
the upper critical dimension dc = 4 the degree of spatial freedom for propagating lines 
is so high that the probability to merge tends to zero. Thus, the contribution of closed 
loops can be neglected and hence the full propagator is effectively described by the free 
propagator. Below the critical dimension, however, loops occur more frequently and begin 
to play a significant role (see Fig. |2|) . 

The loops of the backbone may be associated with the Feynman diagrams of Reggeon 
field theory. As shown in Fig. |2^(a)-(c), the backbone can be decomposed into three 
elementary components. The arrow stands for the free propagato r (|14^ ) while the diagrams 
for branching and merging represent the cubic vertices in Eq. ( |143| ), associated with the 
weights ibr/2. Because of self-destruction A ^ 0, free paths have a finite lifetime, as 
expressed by the bare mass k of the free propagator. Consequently, the paths in a given 
configuration of the backbone have to be weighted by their length. Moreover, each closed 
loop carries a weight — r^/4. 

The negative sign for the weight of closed loops can be explained as follows. The pair 
connectedness function is the sum over all backbone configurations b connecting the sites 
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(a) (b) (c) (d) (e) (f) 

Fig. 24: Feynman diagrams of directed percolation. Left: Elementary components of the backbone: 
(a) free propagator, (b) branching vertex, and (c) merging vertex. Right: One-loop diagrams for (d) 
propagator renormalization and (e,f) vertex renormalization. 

(x',t') and (x, t) weighted by their probabihty Pb' 

c(x',t',x,t) = ^Pfe. (149) 

b 

In order to find an recurrence relation for the probabihty Pb, let us consider directed bond 
percolation on a lattice. Obviously, Pb is the weighted sum over all lattice configurations 
compatible with the backbone configuration b. The backbone itself contributes with a 
factor p"*", where p is the percolation probability and Ub denotes the number of bonds 
occupied by the backbone b. Another factor comes from the bonds outside the back- 
bone. This factor can be expressed as the probability that b is not contained in a larger 
backbone b' . Thus, Pb satisfies the recurrence relation 

Pfe=p"Hl-P~"'' E P"'). (150) 

b',bcb' 

Since b' contains always more loops than b, this relation can be used to expand the pair 



connectedness function ( 149 ) in the number of loops. As can be easily verified, the one-loop 
correction carries a negative sign. More generally, it is possible to show by an inclusion- 
exclusion argument that each closed loop contributes with a negative weight. An analogous 



proof for isotropic percolation is explained in detail in the review article by Essam [108|. 



Thus, apart from the negative weight of closed loops, the backbone may be interpreted 
as a graph consisting of Feynman diagrams. This interpretation is possible because in a 
DP process the field t) represents the local density of particles. It should be noted 
that this is not always true. Moreover, various authors prefer to shift the field by its 
mean field expectation value so that the interpretation as a density is no longer valid. For 
this reason we will continue to use the unshifted fields (j) and (j). 

One-loop approximation 

Slightly below the upper critical dimension one-loop diagrams start to contribute to the full 
propagator while higher-order diagrams are still strongly suppressed. In this regime the 
DP process can be approximated by neglecting higher-order loop diagrams. The resulting 
propagator consists of a sum over n concatenated one-loop diagrams with n running from 
zero to infinity. In momentum space the corresponding expression can be written as 
a simple geometric series. By carrying out the integration one obtains an ultraviolet- 
divergent expression. Hence, in order to regularize the propagator, an upper cutoff 17 
has to be introduced in momentum space. Physically this upper cutoff corresponds to the 
lattice spacing of the DP model. In other words, DP needs a lattice; there is no continuum 
theory of DP. 

In order to approximate the critical exponents, we use Wilson's renormalization group 



scheme [189] which consists of two steps (see Fig. ^5|). At first the theory is coarse-grained 
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Fig. 25: Wilson's renormalization group. Step 1: Scaling transformation in momentum space. Step 2: 
Integration in a momentum shell (shaded region). 



by a scaling transformation x Ax with A < 1, leading to a change of the coefficients in 
the effective action and a dilatation of momentum space (including the cutoff fi). In the 
second step the short-range fluctuations are integrated out in a momentum shell. This can 
be done by evaluating the Feynman diagrams in the range Q < k < A and absorbing the 
resulting contributions in the coefficients. The total change of the coefficients determines 
the RG flow and therefore the critical exponents. 

Let us first consider the scaling transformation. Because of the duality symmetry of 
DP (see Sec. 3^) the action S is invariant under the replacement 

0(x,t) ^ -<^(x,-t), 0(x,t) ^ -</.(x,-t) , (151) 

implying that (j) and (p have exactly the same scaling behavior: 

x^Ax, t^A't, (/>(x, t) ^ A^ 0(Ax, A^t) , (f>{x,t) ^ A^ 4>{Ayi, AH) . (152) 

Under this scaling transformation the effective action (142)-( p^3| ) turns into 

So[(j), 4>] = f d'^X dt 4>ijC, t) (rA^St - DI^X+d+z-2^y2 _ ^^J^X+d+z^"j t) , 

t' D' ft' 

Sint[4>,4>] = \ I d'^xdt rAW+^^ <^(x,t)0(x,t)(0(,t) -<A(,^)) ■ (153) 

r' 

Thus, for an infinitesimal dilatation A = 1 + ^, the four coefficients rescale as 

r' = [l + l{2x + d)]T, 

D' = [l + l{2x + d + z-2)]D, (154) 
k' = [l + l{2x + d + z)]K, 
r' = [l + l{3x + d + z)]T. 

In the second step of Wilson's RG procedure the one-loop diagrams are integrated in a 
momentum shell. The propagator is renormalized by diagram (d) in Fig. ^ 

G^\k,u;r = G,Hk,u;)-^ j^dj,,^, G^{^ + k' ,'^ + u:')Go{^- k' J) , (155) 

where '>' denotes integration in the momentum shell 0, < k < Q/A. This equation can 
be rewritten as 



w2 

k" - D"k'^ + iT"LO = k' - D'k^ + ir'to - —jP , (156) 
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where denotes the integral in Eq. (|155| ). Integrating and expanding the result to 
the lowest order in k and u) yields the series (see Appendix ^ 

Therefore, the coefficients in Eq. ( |156D are renormalized to one-loop order by 

THKd 



t" = t' 



^ - Wrin^D - kY ' ^^^^^ 



K = K 



Finally, we have to renormalize the coupling constant F. Because of the duality symme- 
try ( |119| ) the cubic vertices renormalize identically [see diagrams (e) and (f) in Fig. |2l. 
For the cubic vertices it is sufficient to carry out the integration at k = uj = 0: 

4. Glik,u;)Goi-k,-u;)=T' - ^ . (159) 

Adding the changes of the coefficients under rescaling ( 154D and the subsequent shell 
integration (|158|) - (|159|) we obtain the RG flow equations 



d^D = Z.(2x + . + .-2-^^-|^|^), (160) 



diK = K (^2x + d + z 



Akt{DQ? - K 

Two scaling combinations appear in these equations, namely 

^ 16r(L>1^2 _ ' 4^^pf^2 _ ^) • ^-^"3^^ 

Two of the four parameters r, Z?, k, F can be chosen freelyf^. Here we fix the coefficients 
of the spatial and temporal derivatives, i.e., we require r and D to be invariant under RG 
transformations. Thus the first two equations read 

4-e + 2x- 251 = 0, 2-e + 2x + z-Si = 0, (162) 

where d = 4 — e. The RG flow is then described by two differential equations: 

diK = KU-e + 2x + z-S2) = K (2 + Si - 52) , 
diT = T [A- e + 3x + Z-8S1) = F(e/2-65i) . 



^^On the level of lattice models such as the DK model or the contact process, this freedom corresponds 
to choosing the time scale and a point on the phase transition line. 
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Fig. 26: Linearized RG fiow near the fixed point of directed percolation. To reacli tfie fixed point, tlie 
system has to be on the dashed line, i.e., there is one parameter in the model which has to be tuned to 
criticality. 



At the fixed point (SI, S2), k and T are invariant under RG transformations, i.e., 2 + S^ — 
5*2 = and e/2 — 63^ = 0. Therefore, the fixed point is located at 

Sl = e/12, 5; = 2 + e/12. (164) 



Inserting this solution into Eq. ( |162| ) we obtain two of the three critical exponents, namely 
X = —2 + 7e/12 and z = 2 — e/12. Tlie third exponent can be determined by investigating 
the RG flow in the vicinity of the fixed point. Because of Eq. (|l6lj ), the fixed point values 
for K and F are given by 

^2* _ (2D{2A + e) rrry _iD'^T 2 



where we assumed that VL'^ ~ fi^. Close to the fixed point, the RG flow in Eq. ( |163D 



can be linearized. As shown in Fig. 26, the flow is attractive along the dashed line and 
repulsive elsewhere. To first order in e the corresponding Jacobian is triangular. Hence 
its eigenvalues are given by the diagonal elements 



d^K{2 + Si- 52)L=,.,r=r* = 2 - 6/4 + 0(6^) , 
5rF(e/2-65i)U^,.r^r. = -e + 0(6 



2^ (167) 



The positive eigenvalue corresponds to the repulsive eigenvector (dotted line in Fig. 26). 
Since the parameter k plays the role of the reduced percolation probability p — Pc, this 
eigenvalue is equal to z^J^, rendering the third critical exponent. Because of x = —fi/vi. 
and z = v\\/v± we thus obtain the critical exponents 

/? = 1 - e/6 + 0(e2) , 

= 1/2 + e/16 + 0(e2) , (168) 
M| = 1 + e/i2 + 0(e2) . 
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(a) (b) (c) 

Fig. 27: (a) Schematic drawing of a DP cluster, (b) The same cluster with an absorbing boundary, (c) 
The same cluster in a parabola-shaped geometry. 



A two-loop approximation of these exponents (see Eq. ( |138| )) can be found in Ref. |171]. 
Although the two-loop result is quite accurate in 3-1-1 dimensions, it cannot compete with 
numerical methods in lower dimensions. For example, in 1-1-1 dimensions the approxi- 
mation for the density exponent (3 differs from the known numerical value by more than 
40%. Even one-dimensional fermionic field theories, which have been introduced recently 
in Ref. [19C|, turn out to be inaccurate. Therefore, regarding quantitative results, field- 
theoretic methods are only of limited interest. However, in many cases they are extremely 
useful to understand essential universal properties of the system. For example, various 
scaling relations can only be proven by means of field-theoretic considerations. In fact, the 
field-theoretic renormalization group is one of the most powerful tools of nonequilibrium 
statistical mechanics. 



3.6 Surface critical behavior 

As in equilibrium statistical mechanics, nonequilibrium critical phenomena depend cru- 
cially on the boundary conditions of the system. Because of long-range correlations, the 
choice of the boundary conditions may affect the physical properties of the entire system. 

The critical behavior at surfaces of equilibrium models has been studied extensively (for 
a review see Igloi et al. [|191[| ) . As suggested by Cardy |192(| , surface critical phenomena may 
be described by introducing an additional surface exponent for the order parameter field 
which is generally independent of the other bulk exponents. A similar picture emerges 
in nonequilibrium statistical physics. However, since in this case there is no symmetry 
between space and time, we have to distinguish between spatial, temporal and mixed 
surfaces. The simplest example of a spatial surface is a semi-infinite system with a wall. 
Close to the wall the scaling behavior of the order parameter is characterized by a surface 
critical exponent f5s whose value depends on the type of boundary condition. The most 
important example of a temporal surface is the initial state of a nonequilibrium system. 
As shown below, correlations in the initial state may in fact change the entire evolution 
of a stochastic process. Finally, we will consider systems with mixed boundary conditions 
such as DP in a parabola-shaped space-time geometry. Mixed boundary conditions may 
be viewed as moving boundaries, i.e., the system size varies with time. 
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DP with an absorbing wall 

In DP an absorbing wall may be introduced by cutting all bonds crossing a given (d-l)- 



dimensional hyperplane in space (see Fig. 27). Hence for p > pc the stationary density of 
active sites close to the wall p'f'^^ is expected to be smaller than the density in the bulk. 
In fact, the density at the wall is found to scale as 

Pf"*~(p-Pc)^^ (169) 

with a surface critical exponent Ps > The problem of an absorbing wall was first studied 
in the simpler case of CDP where a surface exponent /3f = 2 was found [193, 194]. In a 



series of papers this scaling theory was later applied to DP with an absorbing wall 1 195 -198] 



(for a review see |199[ ). By means of series expansions and numerical simulations it was 
observed that the mean survival time T of a cluster in the inactive phase next to the wall 
scales as T A^"^^, where A denotes the distance from criticality. In 1+1 dimensions 
the exponent was estimated by 1.0002(3), leading to the remarkable conjecture Ps = 



1 [196]. However, very recent series expansions favor the value 1.00014(2) ^ 1, 



disproving the conjecture ]20C]. In fact, in view of dimensional analysis it seems to be 
unlikely that Ps and i^y are related by a simple linear scaling relation. Moreover, in 
2+1 dimensions the numerical value = 0.26(2) cannot be simply related to the other 
exponents. Similarly, the field-theoretic one-loop result ]198| 

r, = -l/2 + lle/48 + 0(e2), = 3/2 - 7e/48 + 0(6^) (170) 

indicates that the surface exponent is generally independent of the other exponents ([1681). 



198]. A 



The field-theoretic analysis was also extended to DP with an absorbing edge 
closely related application is the study of spreading processes in narrow channels [ 201 ] . It 
is also interesting to study DP with an active wall. This case is related to the problem of 
local persistence and will be discussed below. 

DP clusters in a parabola 

Several years before the problem of an absorbing wall was investigated, Kaiser and Turban 



considered the much more complicated problem of DP in a parabola-shaped geometry P02| , 
203]. Assuming an absorbing boundary of the form x = ^ct" they proposed a general 
scaling theory. It is based on the observation that the width c of the parabola c scales 
as c ^ A^°"~^c under rescaling (p2|). Therefore, the boundary is relevant for a > 1/z and 
irrelevant otherwise. To implement this scaling theory, the scaling forms (Il0l| )- (|102D have 
to be extended by an invariant argument of the form t"~'^l^ jc. The survival probability 
of a cluster (ll02| ), for example, has to be generalized by 

P(i) ~ g{A i^/^ll , f^/'/N, r^-^/Vc) . (171) 

This scaling form is supported by numerical results and a mean field approximation ]l203(] . 
The authors also derived a conjecture for the fractal dimensions 

d||(f7) = 1 - ZCT(d|| - 1) , d±{a) = dii{a)/a , (172) 

where = (/3 + 7)/z^[|. 
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Fig. 28: "Critical initial slip" of the particle density measured in a (1+1) -dimensional directed bond 
percolation process at criticality on a lattice with 10* sites, (a): Particle density p{t) for various initial 
densities po as a function of tim e. T he dashed lines indicate the slopes +9 and —S. (b): Data collapse of 



the same data according to Eq. (175) 



Early-time behavior and critical initial slip 

In Sec. ( |3.4D we reviewed two Monte Carlo techniques for systems with phase transitions 
into absorbing states which differ in their initial state. In simulations starting with a fully 
occupied lattice the particle density at criticality decreases as pit) ~ . On the other 

hand, in dynamic simulations starting from a single particle (active seed), we observe an 
increase of the average number of particles as N(t) ~ t^. In general the exponent 6 is 
independent from the bulk exponents In the case of DP, however, the duality 

symmetry under time reversal (see Eq. (|119D) implies the additional hyperscaling relation 



e = {dv^ - 2(5)/v\\ . (173) 

An interesting crossover phenomenon between initial increase and asymptotic decay of the 
number of particles emerges when a critical spreading process starts with a low-density 



distribution of active sites. Fig. 28 a shows the temporal behavior of the density of active 
sites p{t) for various initial densities po- The density first increases as p{t) ~ until it 
reaches a maximum value at time tc when it crosses over to the usual asymptotic decay 
p(t) ~ This phenomenon is sometimes referred to as the critical initial slip of 

nonequilibrium systems. As can be seen in Fig. ^a, the curves converge to a single one 
after sufficiently long time when the memory of the initial condition is lost. The crossover 
time tc depends on the initial density pQ and scales as 

tc^Po ■ (174) 

In finite-size systems near criticality the critical initial slip may be described by adding 
the scale-invariant argument pQt^^'^^^^^ to the scaling form (101), i.e. 



p{t) ~ t-^/^li /(At^/^li, f^/^N, pot^/^ll+'') • (175) 

The scaling function / behaves asymptotically as /(0, 0, u) ~ n for u ^ and /(0, 0, u) = 
const for u — > oo. To verify this scaling form at criticality, we have plotted p{t)t^ versus 
pot^+^ in Fig. lib. As can be seen, we obtain a convincing data collapse. 



66 



The critical initial slip in a DP process can be interpreted as follows. In a low-density 
initial state the active sites are separated by empty intervals of a certain average size 
^0- As time evolves, they generate individual clusters of connected sites (see Sec. 3.1). 



Initially these clusters are spatially separated; they do not interact and the particle number 
therefore increases as . Only a fraction of these clusters survives, each of them 
spanning a volume of These surviving clusters start touching each other when ~ 

Po^t^c- Therefore, we expect the crossover to take place at tc ~ Pq^^^ '^/^) ^ Insertion of 



the DP hyperscaling relation ( |173| ) leads to Eq. (174). It is worth being mentioned that 
dynamic simulations starting from a single particle represent the limit /?o — > 0. In this 
case tc diverges and the critical initial slip extends to the entire temporal evolution of the 
system. 

Correlated initial conditions 

The previously discussed early-time behavior shows that initial states with short-range 
correlations may affect the temporal evolution of a DP process for a limited time until the 
system crosses over to the usual decay of the particle density. Let us now turn to initial 
states with long-range correlations of the form 

{sis^+r) - r'^-'^ , (176) 

where < cr < d controls the power-law decay of the correlations on large scales. In 
one dimension such states can be generated by creating uncorrelated empty interval^ 
of length i which are algebraically distributed as P{i) ~ There are, however, 

many possibilities to create such states because higher order correlations can be chosen 
freely. Apart from cutoffs, the resulting particle configurations do not exhibit a specific 
length scale .^O) instead they are characterized by a fractal dimension df = a. It turns out 
that long-range correlations may change the entire temporal evolution of a DP process 
(similar phenomena can be observed in other nonequilibrium critical systems such as in 
the annihilation model pO£ 



For a = d the particles are homogeneously distributed, leading to the usual long-time 
behavior p{t) ~ t^^/'^ll . For cr — > the fractal dimension tends to zero, corresponding 
to isolated particles where we expect the density to increase a p{t) t^ . In between 
numerical simulations suggest that the decay exponent changes continuously by 

pW~P„*"<-', „ , (ITT) 

y^yd - a - (i/vx_) for cj > (Tc 

where ctc = (3 /vi. plays a role of a critical threshold above which correlations in the initial 
state become relevant (see Fig. 29). Below cTc the initial distribution of particles is so 



sparse that interactions between growing clusters turn out to be irrelevant. 

The numerical result can be proven by a simple field-theoretic calculation |204|. In 



order to take the initial state into account, the field-theoretic action ( |142[ )-( p!43[) has to be 
extended by the term 



Sic = p I d'^x0(x,O)0o(x) . (178) 



12 



The generation of fractal distribution is crucial since for a < d the p artic le density is zero. Therefore 



appropriate cutoffs have to be introduced, as described in detail in Ref. [204 
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Fig. 29: Correlated initial conditions in a DP proc ess. Numerical estimates for the exponents a{a) in one 
spatial dimension. The theoretical prediction (177) is shown by the solid line. The dashed line indicates 
the 'natural' correlations of DP. 



Here the initial particle distribution is represented by a field (j)o{x) that is coupled to 
the 'creation operator' cj) at time t = via a coupling constant fi. The scaling behavior 
of </>o(x) depends on the fractal dimension. Obviously, homogeneous initial conditions 
00 (x) = const are invariant under rescaling, whereas a fully localized initial condition 
(/)o(x) = 5'^(x) has the scaling dimension —d. Therefore, the initial density should scale as 

M^) - ^"'-"M^) = A"-''</.o(x) . (179) 

It is easy to verify that the contribution Sic does not lead to additional loop corrections in 
the field theory. Therefore, the coupling between the system and the initial condition will 
not be renormalized. Moreover, it can be shown that higher-order contributions of the form 
0'^(x, 0) 0o(x) with k > 1 are irrelevant under renormalization [ 130| , 204 1. Consequently, 
the coupling constant fi scales as 

^ ^ A'^+^fi , (180) 

where x = ~P/'^±- Hence the correlations in the initial state are relevant if cr > fic = 
Scaling invariance of the expression p(t) ~ p^t"^"^ implies that az = d — a + Xi completing 



the proof of Eq. (177). 



Interestingly, the above calculation does not depend on the specific form of correlations 
in the initial state but only on the scaling dimension of the distribution. That is, no matter 



how the particles are distributed - as long as the distribution scales as in Eq. ( 179 ), the 



particle density decays according to the scaling form (177). Moreover, it is interesting to 
note that a critical DP process itself generates two-point correlations (sjSj+r) ~ r~^^'^^, 
corresponding to the 'natural' fractal dimension df^± = d — (3/v±. Therefore, choosing 
'natural' correlations a = d — (3/v± the number of particles remains almost constant (see 
dashed lines in Fig. ^). Similar phenomena have been observed in the Glauber-Ising 
model with correlated initial conditions ||20(: 
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Persistence probability in a DP process 

In the past few years it has been realized that certain first passage quantities of critical 
nonequilibrium processes exhibit a power law decay with non-trivial exponents. One of 
these quantities is the local persistence probability Pi{t), defined as the probability that 
a local variable Si{t) at a given site i does not change its state until time t during the 
temporal evolution. In various systems it was found that 

Piit)^t-'^', (181) 



where 9i is the so-called local persistence exponent p07| - pi3|| . A similar quantity, the global 
persistence probability Pg{t), which is defined as the probability that the global order 
parameter does not change its sign up to time t, is also found to decay as a power law with 



a global persistence exponent Og \21A~2\t\. In general the exponents 6i and Og are different 
and independent of the usual scaling exponents. Since the persistence probabilities depend 
on the history of evolution as a whole|^, it is generally hard to determine these exponents 



analytically. In fact, only a few exact results have been obtained so far 120^, ^ , ^1^ 



Persistence exponents are known to exhibit certain universal properties. For example, the 
local persistence exponent of the two-dimensional Glauber model in the ordered phase 



T < Tc does not depend on T [217-219|, whereas it is non-universal with respect to the 
initial magnetization |210| ]. Most researchers believe that persistence exponents are to 
some extent 'less universal' than ordinary bulk exponents. 

In a DP process the local persistence probability Pi{t) may be defined as the probability 
that a particular site never becomes active up to time t. Numerical simulations suggest 
that the local persistence exponent is given by 

01 = 1.50(1) (182) 

independent of the initial density of active sites |220| . Moreover, the numerical data are 
in agreement with the scaling form 

Pi{t, N, A) ~ t-^' /(Ai^/^'ll , N-^'/H) , (183) 

where A = p—pc denotes the distance from criticality. The local persistence probability can 
also be related to certain return probabilities in a DP process with an absorbing boundary 



or an active source [220|. Similar measurements of the global persistence probability Pg{t) 
suggest that 6g > 9i in agreement with recent field-theoretic results [221|. 



3.7 The influence of quenched disorder 

In DP models it is usually assumed that the percolation probability does not vary in space 
and time. However, in realistic spreading processes the rate for offspring production is 
not homogeneous, rather it fluctuates in space and time. For example, the local density of 
open channels in a porous rock will vary because of inhomogeneities of the material. Sim- 
ilarly, most spreading processes take place in inhomogeneous environments. It is therefore 
important to investigate how quenched disorder affects the critical properties of a spread- 
ing process. It turns out that even weak disorder may affect or even destroy the critical 
behavior of DP. 

^■^Persistence probabilities may be regarded as autocorrelation functions involing infinitely many points. 
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In the DP Langevin equation ( |139| ) the parameter k plays the role of the percolation 
probability. Quenched disorder may be introduced by random variations of k, i.e., by 
adding another Gaussian noise field 77: 

K^K + r]{x,t). (184) 

Thus, the resulting Langevin equation reads 

dtp{x,t) = Kp(x,t)-Ap2(x,t) + L>vV(x,t)+C(x,t)+/9(x,t)7?(x,t). (185) 



The noise ij is quenched in the sense that quantities like the particle density are averaged 
over many independent realizations of the intrinsic noise C while the disorder field r] is 
kept fixed. In the following we distinguish three different types of quenched disorder, 
namely spatially, temporally quenched, and spatio-temporally quenched disorder. The 
three variants of quenched disorder differ in how far they affect the critical behavior of 
DP. 



Spatially quenched disorder 

For spatially quenched disorder, the disorder field r] is defined through the correlations 



r/(x)r/(x') = 7(5 (x - X 



(186) 



where the bar denotes the average over many independent realizations of the disorder field 
(in contrast to the ensemble average (...) over realizations of the intrinsic noise Q. The 
parameter 7 is an amplitude controlling the intensity of disorder. In order to find out 
whether this type of noise affects the critical behavior of DP, let us again consider the 
properties of the Langevin equation under rescaling. At the critical dimension = 4 the 
additional term prj scales as 



pr] 



(187) 



i.e., spatially quenched disorder is a marginal perturbation. Therefore, it may seriously 
affect the critical behavior at the transition. The same result is obtained by considering the 
field-theoretic action. Without quenched noise, DP is described by the action of Reggeon 
field theory (see Sec. |37 



5 



d'^x dt ( 



dt 



(188) 



where (p{x, t) represents the local particle density while (/)(x, t) denotes the Martin-Siggia- 
Rosen response field. As shown by Janssen ]222| ], spatially quenched noise can be taken 
into account by adding the term 



5 + 7 



d'^x 



1 2 



dti 



(189) 



By simple power counting one can prove that this additional term is indeed a marginal 
perturbation. Janssen showed by a field-theoretic analysis that the stable fixed point is 
shifted to an unphysical region, leading to runaway solutions of the flow equations in the 
physical region of interest. Therefore, spatially quenched disorder is expected to crucially 
disturb the critical behavior of DP. The findings are in agreement with earlier numerical 
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Fig. 30: (l+l)-dimensional DP with spatially quenched disorder. Left: In the glassy phase quenched 
disorder forces active sites to percolate in narrow 'channels' where the local percolation probability is high. 
Right: Supercritical disordered DP process starting from a single seed, leading to avalanches (marked by 
the arrows) where the spreading agent overcomes a local barrier. 



results by Moreira and Dickman |223[| who reported non-universal logarithmic behavior 
instead of power laws. Later Cafiero et al. |224| showed that DP with spatially quenched 



randomness can be mapped onto a non-Markovian spreading process with memory, in 
agreement with previous results. 

From a more physical point of view, spatially quenched disorder in 1+1 dimensional 
DP was studied by Webman et al. [ p25| . It turns out that even weak disorder drastically 
modifies the phase diagram. Instead of a single critical point one obtains an intermediate 
phase of very slow glassy-like dynamics. The glassy phase is characterized by non-universal 
exponents which depend on the percolation probability and the disorder amplitude. For 
example, in a supercritical 1+1 dimensional DP process without quenched disorder the 
boundaries of a cluster propagate at constant average velocity v. However, in the glassy 
phase V decays algebraically with time. The corresponding exponent turns out to vary 
continuously with the mean percolation probability. The power-law behavior is due to 



'blockades' at certain sites where the local percolation probability is small (see Fig. 30). 
Similarly, in the subcritical edge of the glassy phase, the spreading agent becomes localized 
at sites with high percolation probability. For d > 1, however, numerical simulations 
indicate that a glassy phase does not exist. 

Temporally quenched disorder 

Temporally quenched disorder is defined by the correlations 



r/(f)r?(t')=7<5(t-0- (190) 

In this case the additional term is a relevant perturbation which scales as pr] — > K~^/'^~^ prj. 
Therefore, we expect the critical behavior and the associated critical exponents to change 
entirely. In the field-theoretic formulation this corresponds to adding a term of the form 



S + -i dt 



n 2 



(191) 



The influence of temporally quenched disorder was investigated in detail in Ref. |llf:]. Em- 



ploying series expansion techniques it was demonstrated that the three exponents /?, u\\ 
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Fig. 31: Spreading process in a coupled map lattice for r = 3 and D — 0.57. Chaotic and laminar sites are 
represented by black and white pixels, respectively. Deterministic synchronous updates lead to symmetric 
clusters (left) with nonuniversal behavior. Asynchronous updates (right) destroy these correlations, leading 
to typical DP clusters. 

vary continuously with the disorder strength. Thus the transition no longer belongs to 
the DP universality class. 

Spatio-temporally quenched disorder 

For spatio-temporally quenched disorder, the noise field rj is uncorrelated in both space 
and time: 

?7(x, t)r/(x', t') = 7 (5'^(x - x')S{t - t') . (192) 
In Reggeon field theory, this would correspond to the addition of the term 

(193) 

being an irrelevant perturbation. In fact, this noise has essentially the same properties as 
the intrinsic noise and can be considered as being annealed. Spatio-temporally quenched 
disorder is expected in systems where each time step takes place in a new spatial environ- 
ment of the system. Examples include water in porous media subjected to a gravitational 
field as well as systems of fiowing sand on an inclined plane (see Sec. ^!g| ). In these cases 
the critical behavior of DP should remain valid on large scales. 



3.8 Related models 

Directed percolation plays a role in various other contexts such as in coupled map lattices, 
the problem of friendly walkers, real-valued spreading processes, models with particle 
conservation, and even in systems with infinitely absorbing states. In the following we 
discuss some of these related models. Moreover, we investigate the special case of compact 
directed percolation in more detail. 

Spreading transitions in deterministic systems 

Spreading transitions can also be observed in certain deterministic lattice models. Instead 
of using random numbers, these models employ chaotic maps in order to generate random 



S ^ S + j / d'^xdt 
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behavior. A simple example of such a chaotic map is given by 



n(t + l) = /(n(t)) , f{x) 



rx 



r(l 



X 



if < X < 1/2 , 
if 1/2 < j; < 1 , 
if 1< X < r/2 , 



(194) 



where r is a free parameter. The chaotic motion of / for x < 1 is governed by a tent map 
of slope r. However, if r exceeds the value 2, the map eventually reaches an absorbing 



state with x > 1, the so-called 'laminar' state of the model. In a coupled map lattice [226] 
many of these local maps Ui{t) are coupled by a diffusive interaction of the form 



Ui{t + 1) = f{u^{t)) + 



D 
'2 



f{u^^lit))-2f{ui{t))+f{u,+l{t)) 



(195) 



where D plays the role of a diffusion constant. The coupled map lattice evolves deter- 
ministically by synchronous updates. By varying D it exhibits a nonequilibrium phase 
transition from a 'chaotic' phase into a 'laminar' state. The existence of absorbing states 
led Pomeau to the conjecture that the transition should belong to the DP universality 



class [227 1, hoping that the apparent randomness of the chaotic maps would effectively 
lead to stochastic spreading of activity on large scales. However, subsequent numerical 
simulations did not agree with this conjecture | 228 |, in particular the exponents were found 
to depend on r. The non-universal behavior of spreading transitions in deterministic sys- 
tems is caused by subtle correlations emerging as artifacts of the deterministic update 
rule. For example, a 'cluster' of chaotic sites starting from a single active seed remains 
symmetric throughout the whole temporal evolution, leading to a qualitatively different 
spreading behavior (see Fig. 31). The consequences of these correlations are not yet fully 
understood. However, replacing the synchronous dynamics of Eq. ( |195| ) by asynchronous 
updates, the deterministic correlations are destroyed and the resulting phase transition 
is indeed characterized by DP exponents |143(| . Similar transitions of two-dimensional 
coupled map lattices have been investigated in Ref. |229]. 



DP and the problem of 'friendly walkers' 

The so-called problem of 'friendly walkers' is defined as follows. Consider the paths of m 
random walkers on a diagonal square lattice. All walks originate in site (0, 0) and end in 
site (x, t). While traveling the walkers may share the same bonds but they are not allowed 
to cross each other (see Fig. |3^ ). In the partition sum Z^i'x, t) each possible configuration 
of random walks is weighted by a factor p'^ , where p > is a parameter and k denotes the 
number of bonds used by at least one of the walkers. For p < 1 it is therefore advantageous 
for the walkers to be 'friendly' to each other, i.e., to share the same bonds. 



Some time ago, Arrowsmith and Essam [P30| suggested a close relationship between 
DP and the problem of friendly walkers. More precisely, they showed that the partition 
function Zm(x, t) is related to the pair-connectedness function c(x,t) of a directed bond 



percolation process by |231] 



c(x, t) = lim Zm{x,t 



(196) 



where p is the usual percolation probability. Here the limit m ^ has to be performed as 
a suitable continuation of polynomial expressions. For example, let us consider m friendly 
random walkers traveling from the origin (0, 0) to the point (x, t) = (0, 2) (see right part of 
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(0,0) 




(x,t)=(0,4) 

Fig. 32: The problem of friendly walkers on a diagonal square lattice. Left: A particular realization of 
three random walkers traveling from the origin (0, 0) to the site (0, 4). Right: All m + 1 realizations of m 
random walkers traveling from (0, 0) to (0, 2). 



Fig. ^). There are m+1 possible configurations; two of them use only two bonds while the 
others use four bonds. Hence the partition function is given by Zm{0, 2) = (m— l)p'^ + 2p^. 
Inserting m = we obtain Zo(0,2) = 2^^ In fact, this expression is exactly equal 

to the pair connectedness function c(0, 2) in a directed bond percolation process. This 
equivalence holds for any (x, t) and also in higher dimensions. Recently this result has 



been generalized to friendly walkers with arbitrary interactions [232|. 



The problem of friendly walkers may also be interpreted as a flow of integer numbers 
on a diagonal square lattice. At the origin there is a source creating an integer number m. 
While traveling on the directed lattice, this integer number may split up into several parts. 
Finally there is a sink where all integers merge into a single one and disappear. Clearly, 
the integers represent just the number of friendly walkers sharing the same bond. 

Even more remarkably, it has been shown that DP is related to the partition sum 
of a chiral Potts model [ p30 , p31 , 233 1, generalizing the well-known result of Fortuin and 
Kasteleyn for isotropic percolation |234]. However, since the definition of the chiral Potts 
model is rather cumbersome, this relation is not of immediate practical benefit. 



DP with real-valued degrees of freedom 

DP models are usually defined in terms of discrete local variables = 0, 1 representing 
inactive and active sites. An interesting variant of DP is 'self-organized directed percola- 
tion', where real- valued local degrees of freedom are used ||235H238(| . To understand the 
basic mechanism, let us consider directed bond percolation. Clearly, a given path between 
two sites is conducting if all bonds along the path are open, i.e., all random numbers 
generated along the path have to be larger than p. Thus, in order to find out whether 
a path is conducting, it is only necessary to keep track of the smallest random number 
generated along this path. This number may be considered as the weight of the path, 
being a measure of its weakest link. However, a pair of sites can be connected by many 
different paths. For the target site to become active, at least one of these paths has to 
be conducting. Therefore, two sites are connected if the maximum of all weights is larger 
than p. 

Interestingly, the maximal weight can be computed by a local update rule which is 
defined in terms of real- valued degrees of freedom Xi{t) £ [0, 1]. Starting with the initial 
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Fig. 33: Self-organized directed percolation in 1+1 dimensions. Left: Spatial configuration of the chain 
after 10* updates. The arrow indicates the percolation threshold of directed bond percolation. Right: 
Minimum value of all sites as a function of time, approaching the percolation threshold of bond DP 
(dashed line). 



condition Xj(0) = the system evolves according to 

Xi{t + 1) = min(^max(Zj~,Xj_i(t)),max(z+,Xj_i(t))^ , (197) 

where we used the notation of Eq. (pOj). A typical spatial configuration of a (1+1)- 
dimensional chain after 10^ updates is shown in Fig |3^. Using this update rule, the binary 
state Si{t) of the corresponding directed bond percolation process can be retrieved by the 
projection 



(198) 



where denotes the heaviside step function. Remarkably, the update rule ( |197| ) does 
not involve the percolation probability p. Instead, it processes all values of p at once 



until a particular value of p is selected by application of the projection rule ( 198 ). Thus, 
'self-organized directed percolation' can be used as a tool for very efficient off-critical 
simulations. 



Spreading process vi^ith particle conservation 

Recently Broeker and Grassberger p39|| introduced another interesting 'self-organized' 
variant of DP which is motivated as follows. A gardener takes care of N plants in a 
flowerbed. The flowers are seized with a parasite. Once a plant is struck, it perishes 
irreversibly. Moreover, the parasite may spread to neighboring plants. However, if the 
number of befallen plants exceeds a certain number, the gardener replaces one of them, 
keeping the number of infected plants constant. In more technical terms, the number of 
active sites is conserved by means of a global update rule. The update consists of two 
steps. At first one of the active sites activates a randomly chosen neighbor, modeling the 
spreading of the parasite. If this move was successful, another randomly chosen active 
site is deactivated, representing the global control of the gardener. Clearly, the number of 
active sites M is conserved, i.e., the model has no absorbing states. 

The number of active sites M is specified by the initial state. For example, we may 
start with a compact domain of M active sites on an infinite lattice. Initially, spreading 
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occurs only at the edges of the domain. As time proceeds, the distribution of active sites 
becomes more and more sparse, forming a diffusing cloud. Nevertheless, the cloud keeps 
its integrity and reaches a typical size after some time. Amazingly, the dynamic processes 
in the interior of the cloud are those of an almost critical DP process. In fact, as shown 
in Ref. ]239| ], most properties of the cloud can be explained in terms of DP scaling laws. 
Considering a small region in the interior of the cloud, the relation to DP is quite obvious: 
The two processes for offspring production A — > 2A and self-destruction A ^ occur 
randomly in space, just as in a contact process with random-sequential updates. However, 
the global control adjusts the ratio of the effective rates and drives the system to criticality. 

Branching Potts interfaces 

Recently Cardy [ p40| studied a field-theory for branching interfaces between ordered 
domains of a g-state Potts model. In two spatial dimensions these interfaces are one- 
dimensional objects. For q < Qc they become fractal with a vanishing interfacial tension 
at the critical point, while for q > Qc the interfacial width diverges at a finite value of 
the tension, indicating a first-order transition. In a certain limit, namely q — > — oo, the 
model becomes equivalent to a DP process. Therefore, the model provides a field theory 
of directed percolation that differs from the standard field theory discussed in Sec. |3.5| . 
Although both field theories 'intersect' in one dimension, they are completely different. 
In particular, the loop expansion starts out from different critical dimensions, namely 
dc = 4 for Reggeon field theory and dc = 2 for branching Potts interfaces. Consequently, 
in the latter case the one-loop estimates for the critical exponents in d = 1 are much more 
accurate. 



DP models with infinitely many absorbing states 

According to the DP conjecture, phase transitions into a single absorbing state belong 
generically to the DP universality class. However, DP behavior may also be observed 
in models with several or even infinitely many absorbing states. An interesting example 
is the dimer-trimer model for heterogeneous catalysis introduced by Kohler and ben- 
Avraham |241]. This model generalizes the ZGB model and is defined by the reaction 
scheme 



00 ^ AA 
000 ^ BBB 
AB 00 



at rate p , 
at rate \ — p 
at rate oo . 



(199) 



On an infinite lattice this model has infinitely many absorbing states. For example, con- 
figurations of dimers and trimers separated by single vacant sites are absorbing. The 
dimer-trimer model displays a phase transition in 2+1 dimensions. Initially, the values 
of the critical exponents were found to be different from those of DP. Later refined sim- 
ulations confirmed, however, that the dimer model still belongs to the DP universality 



class 1 242]. The same result was found in a similar model for catalysis of dimers and 
monomers 1 243, 244]. Another important example is the pair contact process without dif- 



fusion ]245] which is defined by the reaction scheme 

2A^2,A, 2A- 



0. 



(200) 



In this model solitary particles neither react nor diffuse. Starting from random initial 
conditions, the critical pair contact process evolves into certain frozen configurations, as 
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Fig. 34: The pair contact process 2A 3A, 2A ^ at criticality. 



demonstrated in Fig. 34. As can be seen, it is important that single particles are not 
allowed to diffuse. In fact, by adding diffusion of individual particles the critical behavior 



of the model changes entirely (see Sec. 4.5). 

In all models with infinitely many absorbing states and non-conserved order parameter 
the critical exponents /3, i/y , i'± coincide with those of DP. This observation suggests an 
extension of the DP conjecture to systems with several absorbing states which are char- 
acterized by a non-conserved single-component order parameter [246]. However, it was 
realized that the dynamic exponents 6 and 9 depend on the initial condition and even 
violate the usual DP hyperscaling relation ( |120 ). Mendes et al. [129| resolved this prob- 
lem by introducing the generalized hyperscaling relation ( |100| ). However, the sum S + 9 is 
believed to be independent of the initial condition. 

Recently Muhoz et al. proposed a Langevin equation for systems with infinitely many 
absorbing states |147]. It differs from the usual DP Langevin equation (|139| ) by an addi- 
tional term: 



dt(f){x, t) =K(j){x, t) - A</)2(x, t) + L»VV(x, t) + C(x, t) + 



+ a(j){x, t) exp 



t 

-w I dt' (p{-x.,t') 
/o 



(201) 



Here a and w are certain constants (the noise correlations are assumed to be the same as 
in Eq. ( |14C| )). This Langevin equation is non-Markovian, i.e., it has a temporal memory. 
The memory is local since the integral correlates fluctuations at the same position in 
space. From the physical point of view, this memory encodes the local realization of 
the absorbing state. As can be seen in Fig. 34, the emerging inactive domains have a 
highly inhomogeneous structure which can be regarded as a fingerprint of the history 
of the spreading process. A detailed numerical analysis of the Langevin equation ( pOl] ) 
confirmed that the exponenets /3, z^y , and u± do belong to the DP class while 5 and 9 vary 
with the density of the initial state [ 247 1 . In order to explain the apparent nonuniversality 
of the spreading exponents, Grassberger developed a simple toy model that grasps the 
main properties of such spreading processes [ p48[ . In this toy model the spreading rate at 
a given site changes irreversibly at the first encounter with the spreading agent. Although 
the model does not involve multiple absorbing states, it displays similar 'nonuniversal' 
properties. 
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Fig. 35: Snapshot of a critical epidemic process with immunisation grown from a single seed after 200 
time steps. Active and immune sites are represented by bold and thin dots, respectively. An zone of 
activity propagates outwards, leaving a cluster of immune sites behind. The resulting cluster belongs to 
the universality class of isotropic percolation. 



Another important example for systems with infinitely many absorbing states is dam- 
age spreading where two copies of a stochastic system evolve under the same realization 
of thermal noise. The concept of damage spreading will be discussed in detail in Sec. ^. 



Epidemic processes with immunization 

As we have seen in Sec. |3.1| , epidemic models without immunization belong generically 
to the DP universality class. In most cases, however, an infected individual becomes 
increasingly immune after recovery, i.e., the susceptibility for a new infection decreases. 
Cardy and Grassberger showed that epidemic models with immunization are in the same 
universality class as dynamic percolation | 24g| , 250| . It is important to note that dynamic 
percolation differs significantly from directed percolation. For example, let us consider 
a spreading process with immunization in 2+1 dimensions starting from an initial state 
where all sites are equally susceptible for infections. If a single site in the center is infected, 
there is a finite probability that the disease will spread. However, since infected sites 
become increasingly immune, a more or less irregular front of activity moves away from 
the origin, leaving behind a certain cluster of immune sites (see Fig. [35| ). The morphology 
of this cluster depends on the percolation parameter. In the supercritical case there is a 
finite probability that the front moves to infinity, whereas in the subcritical regime the 
process stops after some time. At criticality it turns out that the generated cluster has 
the same asymptotic properties as critical clusters of isotropic percolation [109[ (cf. left 



part of Fig. IC). Thus, dynamic percolation can be used as a tool to generate isotropic 



percolation clusters and should not be confused with directed percolation. In particular, 
the critical exponents turn out to be different in both cases. Interestingly, even a small 
degree of immunization suffices for a (2+l)-dimensional epidemic process to cross over 
from directed to dynamic percolation (together with a shift of the critical point). A 
renormalized field theory of dynamic percolation was studied in |251]. 



Compact directed percolation 

Let us finally come back to compact directed percolation (CDP) [|117| which characterizes 
the critical behavior of the DK model at the upper terminal point of the phase transition 
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line pi = 1/2, P2 = 1 (see Fig. 14). The case p2 = 1 is special because there are two 
symmetric absorbing states, namely the dry state si = . . . = sjy = and the entirely wet 
state si = ... = Siv = l- In contrast to DP, CDP has a global Z2 symmetry 



1 



Si 



Pi 



1 



Pi ■ 



(202) 



Since wet sites cannot spontaneously become dry, compact islands of active sites are 
formed. In 1+1 dimensions CDP is fully equivalent to a zero temperature Ising model with 
Glauber dynamics or the voter model ]252[ |. Expressing the dynamic processes in terms 
of kinks X between wet and dry domains, the kinks perform an annihilating random 
walk X + X — > 0. Therefore, (l+l)-dimensional CDP is exactly solvable |113|. The 
corresponding critical exponents are given by 



6-- 



0, 

1/2, 



/3' 



1, 
= 0, 



2, 

: 1. 



(203) 



It should be noted that these exponents do not comply with the usual DP hyperscaling 



relation (120). However, as pointed out in Ref. [253|, they satisfy the generalized hy- 
perscaling relation ( [1001 ). In fact, as can be verified easily, for CDP the backbone of a 
two-point function (see Sec. |3.5| ) is no longer statistically invariant under time reversal. 



Because of the vanishing exponent P, the CDP transition is discontinuous. In fact, for 
pi < 1/2, p2 = 1 the empty lattice is a stable stationary state while for pi > 1/2 the fully 
occupied lattice is stable. Various spreading models display a crossover from CDP to DP. 
In these models, the rate for the reaction A — > is very small. Therefore, clusters appear 
to be compact on small scales. On larger scales, however, clusters break up into several 
active branches, leading to DP behavior in the asymptotic limit. This type of crossover 



has been studied in detail in Refs. 1 254-257] and may also play a role in experiments of 
flowing granular matter (see Sec. 



3.9 Experimental realizations of directed percolation 

So far we have seen that directed percolation is the generic universality class for nonequilib- 
rium phase transitions into absorbing states. In fact, DP seems to be of similar importance 
as the Ising model in equilibrium statistical mechanics. Despite this success in theoretical 
statistical physics, the critical behavior of DP, especially the values of the critical expo- 
nents, have not yet been confirmed experimentally. The lack of experimental evidence is 
indeed surprising, especially since a large number of possible experimental realizations have 
been suggested in the past. As Grassberger emphasizes in a summary on open problems 



in DP p06| : 



"...there is still no experiment where the critical behavior of DP was seen. This 
is a very strange situation in view of the vast and successive theoretical efforts 
made to understand it. Designing and performing such an experiment has thus 
top priority in my list of open problems. ". 

What might be the reason for the apparent lack of experimental evidence? It seems that 
the basic features of DP, which can easily be implemented on a computer, are quite difficult 
to realize in nature. One of these idealized assumptions is the existence of an absorbing 
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Fig. 36: Left: The reaction rate on Pt(210) as a function of Pco (figure reprinted fr om Ref . [259 ). Right: 
STM image of a catalytic reaction on a Pt(lll) crystal (figure reprinted from Ref. |260|). 



state. In real systems, however, a perfect non-fluctuating state cannot be realized. For 
example, a poisoned catalytic surface is not completely frozen, instead it will always be 
affected by small fluctuations. Although these fluctuations are strongly suppressed, they 
could still be strong enough to 'soften' the transition, making it impossible to quantify 
the critical exponents. 

Another reason might be the influence of quenched disorder due to spatial or temporal 
inhomogeneities. In most experiments frozen randomness is expected to play a signiflcant 
role. For example, a real catalytic surface is not fully homogeneous but characterized by 



certain defects leading to spatially quenched disorder. As has been shown in Sec. 3.7, this 
type of disorder may affect or even destroy the critical behavior of DP. 

In the following we summarize some of the most important experimental applications 



which have been proposed so far |255]. Other experimental applications in systems of 



growing interfaces will be discussed in Sec. 
Catalytic reactions 

It is well known that under speciflc conditions certain catalytic reactions mimic the mi- 
croscopic rules of DP models. For example, as shown in Fig. 0, the ZGB model for the 
catalytic reaction CO + O ^ CO2 on a platinum surface displays a continuous transition 
at y = ui belonging to DP. In real catalytic reactions, however, only the discontinuous 



transition at y = 1/2 can be observed. Fig. 36 shows the reaction rates as functions of 



the CO pressure measured in a catalytic reaction on a Pt(210) surface |25gf| . Although 
this experiment was designed in order to investigate the technologically interesting regime 
of high activity close to the flrst-order phase transition, the results clearly indicate that 
poisoning with oxygen does not occur. Instead the reactivity increases almost linearly 
with the CO pressure. Similar results were obtained for Pt(lll) and for other catalytic 
materials. Thus, so far there is no experimental evidence for DP transitions in catalytic 
reactions. 

One may speculate why the DP transition is obscured or even destroyed under exper- 
imental conditions. On the one hand, the reaction chain in the experiment is much more 
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Fig. 37: Simplified drawing of the Douady-Daerr experiment. Pouring sand on top of a plane with 
inclination (j), a thin layer settles and remains immobile. Perturbing the layer locally with a stick we can 
release an avalanche of flowing. 



complicated than in the ZGB model |261]. Moreover, the 0-poisoned system might not 
be a perfect absorbing state, i.e., the surface can still adsorb CO molecules although it is 
already saturated. Another possibility is thermal (nonreactive) desorption of oxygen, act- 
ing as an external field which drives the system away from criticality [ p62| [cf. Eq. ( |116D ]. 
Finally, defects and inhomogeneities of the catalytic material could lead to an effective 
(spatially quenched) disorder. 

For a long time the microscopic dynamics processes were difficult to study experimen- 
tally. However, in recent years novel techniques such as scanning tunneling microscopy 
(STM) led to an enormous progress in the understanding of catalytic reactions, pointing 
at various unexpected subtleties. For example, it was observed that reactions preferably 
take place at the perimeter of oxygen islands ||260(| . Furthermore, it was realized that 
adsorbed CO molecules on Pt(lll) may form three different rotational patterns represent- 
ing the c(4x2) structure of CO on platinum, i.e., there are several competing absorbing 
states p63|] . Moreover, the STM technique allows one to trace individual molecular re- 
actions and to determine the corresponding reaction rates. In addition, the influence of 



defects such as terraces on catalytic reactions can be quantified experimentally |264]. We 
may therefore expect a considerable progress in the understanding of catalytic reactions 
in near future. 

Flowing granular matter 

Recently it has been shown that simple systems of flowing sand on an inclined plane, such 
as the experiments performed by Douady and Daerr [265, p66 |, could serve as experimental 



realizations of DP [ P67| . In the Douady-Daerr experiment glass beads with a diameter of 
250-425 /xm are poured uniformly at the top of an inclined plane covered by a rough velvet 
cloth (see Fig. |37| ). As the beads flow down, a thin layer settles and remains immobile. 
Increasing the angle of inclination cj) by the layer becomes dynamically unstable, i.e., 
by locally perturbing the system at the top of the plane an avalanche of flowing granular 
matter will be released. In the experiment these avalanches have the shape of a fairly 
regular triangle with an opening angle 9. As the increment A99 decreases, the value of 9 
decreases, vanishing as 

tan0~(Av9)^ (204) 
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Fig. 38: Typical clusters generated at criticality on small and large scales, illustrating the crossover from 
CDP to DP. 



with a certain critical exponent x. The experimental results suggest the value x = 1 |266|. 



In order to explain the experimentally observed triangular form of the avalanches, 
Bouchaud et al. proposed a mean- field theory based on deterministic equations, taking 



the actual local thickness of the flowing avalanche into account |268|. This theory predicts 



the exponent x = 1/2. Another explanation assumes that flowing sand may be interpreted 



as a nearest-neighbor spreading process |267|. Here the avalanche is considered as a cluster 



of active sites. Identifying the vertical coordinate of the plane with time and the increment 
of inclination Aip with p — pc, the opening angle is expected to scale as 

tan^ ~ ^x/^ll ~ {Aippi-"^ , (205) 

where i^y and iy± are the scaling exponents of the spreading process under consideration. 
To support this scaling argument, a simple lattice model was introduced which mimics 



the physics of flowing sand [267|. The model exhibits a transition from an inactive to an 
active phase with avalanches whose compact shapes reproduce the experimental observa- 
tions. On laboratory scales the model predicts a transition belonging to the universality 
class of compact directed percolation [see Eq. ( |203| )], implying that 

x = u\\-iy± = 1. (206) 

The CDP behavior, however, is only transient and crosses over to DP after a very long 
time. Thus the Douady-Daerr experiment - performed on sufficiently large scales - may 
serve as a physical realization of DP. Irregularities of the layers thickness can be considered 
as spatio-temporally quenched disorder which is irrelevant on large scales (see Sec. 3.7). 



Thus, in contrast to catalytic reactions, the problem of quenched disorder does not play a 
major role in this type of experiments. 

The crossover from CDP to DP is very slow and presently not accessible in experiments. 



To illustrate the crossover, two avalanches are plotted on different scales in Fig. 38. The 
left one represents a typical avalanche within the first few thousand time steps. As can be 
seen, the cluster appears to be compact. However, as shown in the right panel of the figure, 
the cluster breaks up into several branches after a very long time. Recent experimental 
studies |26£] confirm that for high angles of inclination critical avalanches do split up 



into several branches (see Fig. |39|). Yet here the avalanches have no well defined front, the 
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Fig. 39: Left: Triangular avalanches in the Douady-Daerr experiment. Right: Avalanches splitting up 
into several branches for high angles of inclination (reprinted with kind permission from A. Daerr). 



propagation velocity of separate branches rather depends on their thickness. It is therefore 
no longer possible to interpret the vertical time coordinate. Another problem 

is the kinetic energy of the grains. According to arguments by Dickman et al. p38[ |, 
continuous phase transitions into absorbing states can only be observed if the inertia of 
particles can be neglected. 

Finally, it is not yet known how the spreading process depends on correlations in the 
initial state. As shown in Sec. 3^, such long-range correlations may change the values of 
certain dynamic critical exponents. However, recent studies of a single rolling grain on an 
inclined rough plane |270| support that there are presumably no long-range correlations 
due to a 'memory' of rolling grains. By means of molecular dynamics simulations it was 
shown that the motion of a rolling grain consists of many small bounces on each grain 
of the supporting layer. Therefore, the rolling grain quickly dissipates almost all of the 
energy gain from the previous step and thus forgets its history very fast. For this reason 
it seems to be unlikely that quenched disorder of the prepared layer involves long-range 
correlations. Therefore, flowing granular matter seems to be a promising candidate for an 
experimental realization of DP. 



Porous media 

DP is often motivated as a model for water percolating through a porous medium in a 
gravitational field. Due to an external driving force, the flow in the medium is assumed 
to be strictly unidirectional, i.e., the water can only flow downwards (in contrast to the 
depinning models of Sec. where the water can flow forth and back). Although this 
application seems to be quite natural, it is difficult to realize experimentally. As shown in 
Ref. [ ^7l| , porous media in nature are highly irregular. By cutting sandstone into slices 
and digitizing the section images, the porosity distribution and the local connectivity were 
measured and averaged over 99 samples. As expected, the pores have different sizes and 
are distributed irregularly. In addition, the percolation probability is found to depend on 
the local porosity and the direction in space, i.e., sandstone is an anisotropic material. 
But there are even more fundamental problems. On the one hand, water is a conserved 
quantity, leading to unpredictable long-range correlations in the bulk. On the other hand, 
water can always flow against the gravity field by means of capillary forces. Therefore, 
it is quite difficult or even impossible to verify scaling laws in such experiments and is 
not yet clear whether the relation to DP is meaningful or simply a commonly accepted 
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misconcept. 



Epidemics 

Another frequently quoted application of DP is the spreading of epidemics without immu- 



nization 1 101, 111]. In an epidemic process infection and recovery resemble the reaction- 
diffusion scheme of DP If the rate of infection is very low, the infectious disease will 
disappear after some time. If infections occur more frequently, the disease may spread 
and survive for a very long time. However, spreading processes in nature are usually 
not homogeneous enough to reproduce the critical behavior of DP. Moreover, in many 
realistic spreading processes short-range interactions are no longer appropriate. This sit- 
uation emerges, for example, when an infectious disease is transported by insects. Such 
long-range interactions may be described by Levy flights, leading to continuously varying 



critical exponents (see Sec. |4.lD . 



Forest fires 

A closely related problem is the spreading of forest fires |102, 134[ . Tephany et al. studied 



the propagation of flame fronts on a random lattice both under quiescent conditions and in 
a wind tunnel |272]. The experimental estimates of the critical exponents at the spreading 



transition are in rough agreement with the predictions of isotropic and directed percolation, 
respectively. However, the accuracy of these experiments remains limited. 

Calcium dynamics 

DP transitions may also occur in certain kinetic models for the dynamics of Calcium ions in 
living cells. Ca^^ ions play an important physiological role as second messenger for various 



purposes ranging from hormonal release to the activation of egg cells by fertilization |273, 



274]. The cell uses nonlinear propagation of increasing intracellular Ca^"*" concentration. 



so-called calcium waves, as a tool to transmit signals over distances that are much longer 
than the diffusion length. For example, propagating Ca^+ waves can be observed in 
the immature Xenopus laevis oocyte |275(] . So far theoretical work focused mainly on 
deterministic reaction-diffusion equations in the continuum, explaining various phenomena 
such as solitary and spiral waves ] 276[ . Recently improved models have been introduced 



which take also the stochastic nature of Calcium release into account ]277]. As expected. 



the transition in one of these models belongs to the DP universality class ]^78 ]. However 



from the experimental point of view it seems to be impossible to confirm or disprove this 
conjecture. On the one hand, the size of a living cell is only a few order of magnitude larger 
than the diffusion length, leading to strong finite-size effects in the experiment. On the 
other hand, inhomogeneities as well as internal structures of the cell lead to a completely 
unpredictable form of quenched noise. Therefore, it seems to be impossible to identify the 
universality class of the transition in such experiments. It would be rather an achievement 
to find clear evidence for the very existence of a phase transition between survival and 
extinction of propagating calcium waves. 

Directed polymers 



DP is also related to the problem of directed polymers ]27£]. In contrast to DP, which is 
defined as a local process, the directed polymer problem selects directed paths in a random 
medium by global optimization. Under certain conditions, namely a bimodal distribution of 
random numbers, both problems were shown to be closely related ]|280|. More specifically. 
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the roughness exponent of the optimal path in a directed polymer problem is predicted to 
cross over from the KPZ value 2/3 to the DP value v\\/v± — 0.63 at the transition point. 
Directed polymers were used to describe the propagation of cracks [ p81| . However, in such 
experiments it is usually impossible to verify the tiny crossover from KPZ to DP. 

Turbulence 

Finally, DP has also been considered as a toy model for turbulence. As suggested in 
Ref. [ |227| ], the front between turbulent and laminar flow should exhibit the critical be- 
havior of DP. For example, the velocity of the front should scale algebraically with a 
combination of DP exponents. However, these predictions are based rather on heuristic 
arguments than on rigorous results. In fact, in many respects turbulent phenomena show 
a much richer behavior than DP. Nevertheless there are certain similarities between DP 
and turbulence. Therefore, the study of DP could be helpful for a better understanding 
of turbulent phenomena. 

Summary and outlook 

Directed percolation is keeping theoretical physicists fascinated since more than four 
decades. Several reasons make directed percolation so appealing. First of all, DP is a 
very simple model in terms of its dynamic rules. Nevertheless, the DP phase transition 
turns out to be highly nontrivial. In fact, DP belongs to the very few critical phenomena 
which have not yet been solved exactly in one spatial dimension. Therefore, the critical 
exponents are not yet known analytically. High-precision estimates indicate that they 
might be given rather by irrational than by simple fractional values. Moreover, DP is 
extremely robust. It stands for a whole universality class of phase transitions from a 
fluctuating phase into absorbing states. In fact, a large variety of models displays phase 
transitions belonging to the DP universality class. Thus, on the theoretical level, DP plays 
the role of a standard universality class similar to the Ising model in equilibrium statistical 
mechanics. 

In spite of its simplicity, no experiment is known which confirms the values of the 
critical exponents quantitatively. An exception may be the wetting experiment performed 
by Buldyrev et al. (see below in Sec. |6.2| ), where the value of the roughness exponent a 
coincides with v^/vw within less than 10%. However, since the results of similar experi- 
ments are scattered over a wide range, further experimental effort in this direction would 
be needed in order to confirm the existence of DP in this type of systems. 

Apart from difficulties to realize a non-fluctuating absorbing state, a fundamental prob- 
lem of DP experiments is the emergence of quenched disorder caused by inhomogeneities 
of the system. Depending on the type of disorder, even weak inhomogeneities might ob- 
scure or even destroy the DP transition. Therefore, the most promising experiments are 
those where quenched disorder is irrelevant on large scales. This is the case, for example, 
in wetting experiments and systems of flowing granular matter. 

Although there is certainly a lack of experimental evidence, there is no reason to 
believe that DP is a purely artificial model. To be optimistic, it is helpful to recall the 
history of the Ising model, which has been introduced almost one century ago. Although 
the Ising model is probably the best studied system in equilibrium statistical mechanics, 
there are only few experiments in which the critical exponents have been reproduced (for 
a review see Ref. @). For this reason, many physicists believe that DP should have a 
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counterpart in reality as well, mostly because of its simplicity and robustness. In this 
respect, Grassbergers message remains valid: The experimental realization of DP is an 
outstanding problem of top priority. 
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4 Other classes of spreading transitions 



This Section discusses various other types of nonequihbrium phase transitions into ab- 
sorbing states which do not belong to the universahty class of directed percolation. In 
particular we will address spreading processes with long-range interactions and additional 
symmetries. 



4.1 Long-range spreading processes 

According to the DP conjecture (see Section |3.3D phase transitions in spreading models 
with short range interactions generically belong to the DP universality class. In many 
realistic spreading processes, however, short-range interactions do not appropriately de- 
scribe the underlying transport mechanism. This situation emerges, for example, if an 
infectious disease is transported by insects. Typically the motion of the insects is not a 
random walk, one rather observes occasional flights over long distances before the next 
infection occurs. Similar phenomena are expected when the spreading agent is subjected 
to a turbulent flow. Intuitively it is clear that occasional spreading over long distances will 
significantly alter the spreading properties. On a theoretical level such a super- diffusive 



transport may well be described by Levy flights [104|, i.e., by uncorrelated random moves 



over long distances r which are algebraically distributed as 

P(r) ~ 1//+'^ , (a > 0) . (207) 

The exponent o" is a free parameter that controls the characteristic shape of the distribu- 
tion. The algebraic tale leads to occasional long-distance flights, as shown in Fig. EG. 



Anomalous directed percolation, as originally proposed by Mollison |101| in the context 
of epidemic spreading, is a generalization of DP in which the spreading agent is transported 
by Levy flights. As in the case of ordinary DP, we expect anomalous DP to be characterized 
by certain universal critical exponents /?, z^^, and The question is how these exponents 
depend on a, whether they are independent from one another, and how they cross over to 



the exponents of ordinary DP. Based on field-theoretic considerations, Grassberger | 282 | 
claimed that the critical exponents of anomalous DP should depend continuously on the 
control exponent a. Very recently this work has been considerably clarified and extended 



by Janssen et al. [283|, who presented a comprehensive field-theoretic analysis of anomalous 



spreading processes with and without immunization. 

Anomalous directed percolation: Field-theoretic predictions 

In order to include long-range spreading in the Langevin equation ( |13S| ), the Laplacian 
has to be replaced by a non-local expression. This term can be written as an integral that 
describes Levy flights over the distance r according to the probability distribution P(r): 

at0(x, t) = K (/)(X, t)-\ </.2(x, t) + C(X, t) 

+ D j (fx' P(|x - x'l) [0(x', t) - ^{x, t)] . 

The two contributions in the integrand describe gain and loss processes, respectively. 
Keeping the most relevant terms in a small momentum expansion |283(| , this equation may 
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be written as 

9t0(x, t) = ( DnV + DaVa + ^ ) '/'(x, t) - A(^'(x, t) + C(x, t) , (209) 



where the noise correlations are assumed to be the same as in Eq. ( |140D . D]\[ and Da 
are the rates for normal and anomalous diffusion, respectively. The anomalous diffusion 
operator describes moves over long distances and is defined through its action in 
momentum space 

V^e^'^"' = -Fe^*^-^ , (210) 

where k = |k|. The standard diffusive term D^V"^ takes the short-range component of 
the Levy distribution into account. Note that even if this term were not initially included, 
it would still be generated under renormalization of the theory. The mean-field theory 
of anomalous DP is completely analogous to that of ordinary DP. For a < 2 a scaling 
analysis yields the mean field results 

4 = 2a, /?^^ = 1, vf^ = l/a, i^|f^ = l- (211) 

For (7 > 2 these exponents cross over smoothly to the ordinary DP mean-field expo- 
nents ( |124| ). The mean field approximation is expected to be quantitatively accurate 
above the upper critical dimension dc, while for d < dc fiuctuation effects have to be taken 
into account. By using standard techniques one can derive the effective action 



4>] = j d'^x dt [4>{dt -K- DnV^ - DaV''a)(I) + ^ 



This expression differs from the usual action of Reggeon field theory ( |142| )-( p^ ) by the 
addition of a term representing anomalous diffusion. Simple power counting on this action 
confirms that the upper critical dimension is dc = 2a, below which fluctuation effects 
become important. The field-theoretic RG calculation basically follows the same lines 
as in the case of DP. The resulting critical exponents to one-loop order in d = 2a — e 
dimensions are given by 

f3 = l- 2e/7a + O (e^) , 
i/_L = l/o- + 2e/7a^ + O (e^) , (212) 
= 1 + e/7a + O (e^) . 
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Fig. 41: Critical anomalous directed percolation in 1+1 dimensions for different values of a. The figure 
shows typical clusters starting in the center of the lattice. The case a = 00 corresponds to ordinary DP. As 
a decreases, spatial structures become more and more smeared out until in the mean field limit a = 1/2 
the particles appear to be randomly distributed over the whole system. For small values of a the system 
quickly reaches the absorbing state due to extremely strong finite size effects. 



Moreover, it can be shown that the hyperscahng relation ( |120| ) holds for arbitrary values 
of a. Thus, to one-loop order, 9 and S are given by 

= e/7a + O{e'^) , 5 = 1 - 3e/7(T + O (e^) . (213) 

Finally, since Dj^ will not be renormalized, one can prove the additional exact scaling 
relation 



^11 



iy±{a -d) -2p = . (214) 



This equation implies that anomalous DP is described by two rather than three indepen- 
dent critical exponents. Moreover, it has another surprising consequence. Assuming that 
/3, i'± and z^|| change continuously with a and cross over smoothly, it predicts for fixed d 
the value ctc where the system should cross over to ordinary DP. In order to compute cJc we 
simply have to insert the numerically known values of the DP exponents into Eq. ( pl4|) . 
Surprisingly one obtains cxc = 2.0766(2) in one, Cc — 2.2 in two, and cTc = 2 + e/12 in 
d = 4 — e spatial dimensions. Thus, the crossover takes place at cTc > 2 which collides with 
the intuitive argument that the anomalous diffusion operator should only be relevant 



if cj < 2. But, as pointed out in Ref. [283|, this naive argument may be wrong in an 



interacting theory where the critical behavior is determined by a nontrivial fixed point of 
a RG transformation. The field-theoretic calculation rather predicts anomalous diffusion 
to be relevant in the range 2 < a < add) for d < 4. The surprising conclusion would be 
that 

V\ / (215) 

in certain interacting theories. Loosely speaking, the tendency to correlate particles in 
local spots of activity makes a DP process more sensitive to long-range flights. Therefore, 
the relevance of Levy flights sets in earlier than in the case of simple diffusion. 



A lattice model for anomalous directed percolation 

The field-theoretic predictions can be verified numerically by studying a lattice model for 
anomalous DP that generalizes directed bond percolation |284]. The model is defined on 
a tilted square lattice and evolves by parallel updates. As usual, a binary variable Si{t) 
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Fig. 42: Numerical estimates for the dynamic critical exponents /3, and the derived exponents and v\\ 
in comparison with the field-theoretic predictions (solid lines) and the DP exponents (dot-dashed lines). 

is attached to each lattice site i. Si = \ means that the site is active (infected) whereas 
Sj = denotes an inactive (healthy) site. The dynamic rules depend on two parameters, 
namely the control exponent cr > and the bond probability < p < 1. For a given 
configuration {sj(t)} at time t, the next configuration {si{t + 1)} is constructed as follows. 
First the new configuration is initialized by setting Si{t + 1) := 0. Then a loop over all 
active sites i in the previous configuration is executed. In the (l-l-l)-dimensional case this 
loop consists of the following steps: 




1. Generate two random numbers zl and from a fiat distribution between and 1. 

2. Define two real-valued spreading distances = z^^^^'^ and r/j = z^^" , for spreading 
to the left (L) and to the right (R). The corresponding integer spreading distances 
dL and are defined as the largest integer numbers that are smaller than and 
r/j, respectively. 

3. Generate two further random numbers ul and yn drawn from a flat distribution 
between and 1, and assign Sj_|_i_2d^(t -|- 1) := 1 if < p and Sj_i+2d^(i -|- 1) := 1 
if yR < p. In finite systems the arithmetic operations in the indices are carried out 
modulo N by assuming periodic boundary conditions, i.e. Si = Si±i\f. 



This model includes two special cases. For cr ^ cxd it reduces to ordinary directed bond 



percolation (see Sec. 3.1). On the other hand, for a — > the interaction becomes totally 
random. In this case the mean-field approximation becomes exact with a transition taking 
place at pc = 1/2. In between, the spreading properties of the model change drastically, as 
demonstrated in Fig. As can be verified easily, the assignment r = z~^^'^ reproduces 
the normalized probability distribution 

C (j/j.i+0- jf r > 1 
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Fig. 43: Anomalous annihilation process. The graph shows direct and extrapolated estimates for the 
decay exponent a, as a function of a. The solid line represents the exact result (neglecting logarithmic 
corrections at cr = 1). 



As usual, the distribution has a lower cutoff at r^m = representing the lattice spacing. 



Yet, in contrast to other models |l285| , P86| , no upper cutoff is introduced. In order to reduce 
finite-size effects, the target site is determined by assuming periodic boundary conditions, 
i.e., the particle may 'revolve' several times around the system. 

An interesting aspect of anomalous DP is the possibility to choose a in such a way 
that the critical dimension dc = 2a approaches the actual physical dimension where the 
simulations are performed. Even in one spatial dimension this allows the one-loop re- 



sults ( 212 ) to be verified. For example, if o" = 1/2 + ^, the critical dimension of the system 
is dc = 1 -|- 2/i. Hence the exponents in a (l-l-l)-dimensional system change to first order 
in fi as 

P = 1 - 8///7 + O (^2) ^ ^ ^ l/2 + 5fi/7 + 0{^i^) , 

i^x = 2 - 12///7 + O , S = 1 - 12/i/7 + (/i^) , (217) 

i^ll = 1 + 4^/7 + O , e = 4///7 + 0(/i2) . 

In Fig. the predicted initial slopes are indicated by solid lines. Clearly they are in fair 
agreement with the numerical estimates, confirming the field-theoretic results of Eq. ( pl2|) . 
This is one of the rare cases where we can directly 'see' the field-theoretic results in the 
simulation data. 

Anomalous annihilation process 

The more simple case of anomalous pair annihilation A + A ^ (}> with long-range hop- 
ping [287 1 can be solved exactly by a similar field-theoretic analysis. In the ordinary an- 



nihilation process with short-range interactions, the average particle density is known 
to decay as in Eq. (p7|). The Levy-flight case may be described theoretically by inserting 



an additional operator into the field-theoretic action for pair annihilation [44|. The 



resulting action, which can also be derived systematically |284], reads 



S[(p,4>] = j d'^xdti^4>{dt-DNV^ -DaV''j^)(P + 2X4>(P'^ + X4>'^(I)'^ -(Po4>6{t)^ , (218) 
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where cpo represents the mitial (homogeneous) density at t = 0. An analysis of this action 
follows very closely that of Ref. |4^. For a < 2, power counting reveals the upper critical 
dimension of the model to he dc = cr < 2. For d > dc mean-field theory is expected to 
be quantitatively accurate, with an asymptotic density decay ~ t~^. Below dc, however, 
the renormalized reaction rate flows to an order e = a — d fixed point. The decay of the 
density can therefore be predicted via dimensional arguments (see Fig. 

( t-'^/" for d < fj , 
p{t) ~ <^ Int ioT d = dc = a , (219) 
[ t^^ ioi d> a . 



4.2 Absorbing-state transitions in systems with additional symmetries 



As outlined previously, directed percolation is the canonical universality class for phase 
transitions into a single absorbing states. According to the DP-conjecture, we may there- 
fore expect non-DP behavior to occur in systems with several absorbing states. However, 
it is important to note that the existence of several absorbing states alone does not auto- 



matically lead to non-DP behavior at the transition point. As we have seen in Sec. 3.8, 
even models with infinitely many absorbing states may still belong to the DP universal- 
ity class. Non-DP critical behavior emerges only if there is a symmetry among different 
absorbing states. 

The first examples of non-DP transitions into absorbing states were discovered by 



Grassberger et al. |28^ , 289 1 who observed "a new type of kinetic critical phenomenon" in 
certain one-dimensional stochastic cellular automata. The density exponent (3 ~ 0.6(2) in 
1+1 dimensions was found to differ significantly from the usual DP exponent (3 ~ 0.277. 
Partially because of the complicated dynamic rules of these models it took almost ten 
years until the mechanism behind this type of non-DP behavior was clearly identified. 

Up to now two universality classes of spreading transitions with non-DP behavior 
have been found, namely the so-called parity- conserving class (PC) and Z2-symmetric 
directed percolation (DP2). The PC class is represented most prominently by branching- 
annihilating random walks with even number of offspring (BAWE) [|290| , where the number 
of particles is preserved modulo 2. The DP2 class, on the other hand, which is also referred 
to as the directed Ising class, introduces two symmetric absorbing states. As we will 
see below, both universality classes coincide in one spatial dimension wherefore they are 
usually considered to be identical. However, it is important to note that they differ from 
each other in higher dimensions. 

One may speculate whether systems with a symmetry among several spreading agents 
will also be able to display novel critical properties. However, such multi-color spreading 



processes were found to belong to the DP class as well [291|. A field-theoretic analysis 
confirms this observation [|292|| and predicts that the symmetry among the spreading agents 
may be spontaneously broken. 



The parity-conserving universality class 

The parity-conserving universality class is represented most prominently by branching 
annihilating walks with an even number of offspring |290| , p93| -297|. These non-conserved 
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Fig. 44: The parity-conserving universality class. Typical space-time trajectories of a branching annihi- 
lating random walk with two ofltspring below, at and above the critical point. 

random walks are defined by the reaction-diffusion scheme 

A^{n + 1)A, 2A^ 0, (220) 

A a 

where n = 2, 4, 6, . . . denotes the number of offspring. As an essential feature, this process 
conserves the number of particles modulo 2. Early numerical studies in 1-1-1 dimensions 
assumed instantaneous on-site annihilation a = oo. In this case the model displays a 
continuous phase transition only for n > 4 1 293 -295], while there is no such transition for 



n = 2. Later Zhong and ben-Avraham demonstrated that a phase transition also emerges 
in the case of two offspring, provided that the annihilation rate a is finite |290|. Fig. ^ 



shows a typical cluster grown from a single seed. In contrast to DP, a PC process starting 
with an odd number of particles cannot reach the empty state since at least one particle is 
left. Therefore, initial states with even and odd number of particles are expected to lead 
to different cluster morphologies. 

The relaxational properties of PC models in the subcritical phase differ significantly 
from the standard DP behavior. While the particle density in DP models decays exponen- 
tially as p{t) ~ e~*^''il , in PC models it decays algebraically in the long-time limit. More 
precisely, the temporal evolution of PC processes in the inactive phase is governed by the 
annihilation process 2 A — > 0. Under RG transformations we therefore expect the system 
to fiow towards the fixed point of particle annihilation. Consequently, in the subcritical 
phase the particle density decays algebraically as in Eq. d37 



A systematic field theory for PC models has been presented in Refs. | 296| , 297], con 



firming the existence of such an annihilation fixed point. However, the field-theoretic 
treatment of the (l-l-l)-dimensional case poses considerable difficulties. They stem from 
the presence of two critical dimensions: dc = 2, above which mean-field theory applies, 
and d'^ « 4/3, where for d > d'^ {d < d'^) the branching process is relevant (irrelevant) at 
the annihilation fixed point. Therefore, the physically interesting spatial dimension d = \ 
cannot be accessed by a controlled e-expansion down from upper critical dimension dc = 2. 
Nevertheless the usual scaling theory is still valid below d'^. Currently the best numerical 
estimates of the critical exponents in 1-1-1 dimensions are 

/? = 0.92(2), 5 + 9 = 0.286(2), 

z/||= 3.22(6), 2/z = 1.15(1), (221) 

z/_L = 1.83(3) . 

The actual values of 6 and 9 in dynamic simulations depend on the initial condition. If 
the process starts with a single particle, it will never stop, hence 6 = 0. On the other 
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hand, if the initial seed consists of two particles, one observers that the roles of 6 and 9 
are exchanged, i.e. 6 = 0. 

It has been customary to investigate whether the numerical estimates of the critical 
exponents can be fitted by simple rational numbers. In fact, the estimates of Eq. ( ^2l\ ) 
are in good agreement with the rational values |295[ /3 = 12/13, i^y = 42/13, i'± = 24/13, 
6 + 6 = 2/7, and z = 8/7. In particular, (3/u_i_, the exponent for the decay of spatial 
correlations at criticality, should be equal to 1/2. In the past, however, rational values 
were also proposed for the DP exponents and later disproved by more accurate numerical 
estimates |29ql. 



The DP2 universality class 

The DP2 universality class describes phase transitions in spreading models with two sym- 
metric absorbing states. Since the two absorbing states compete one another, the resulting 
critical behavior is different from ordinary directed percolation]^. In various models, for 



example in certain cellular automata 



as well as in interacting monomer-dimer 



models |29S-302], the two absorbing states emerge as checkerboard-like configurations of 
particles at even or odd sites, respectively (see Fig. ^). Other DP2 models explicitly intro- 
duce two symmetric inactive states, as, for example, nonequilibrium Ising models | 303| , p04 |, 
Z2-symmetric generalizations of the DK model and the contact process [p63|| , monomer- 
monomer surface reaction models with two absorbing states [ |305| ], and certain cellular 
automata |306]. Even in monomer- monomer models with three absorbing states a DP2 
transition emerges at certain points in the parameter space [307,305]. 

On a phenomenological level, spreading transitions with several absorbing states may 
be defined by introducing a single active state A and n symmetric absorbing states 
Ji, . . . , which can be regarded as having different colors. The system evolves according 
to the following descriptive rules: 



1. Spreading of activity: Active sites turn their inactive nearest neighbors into the 
active state. 

^''Without symmetry, one of the two absorbing states will dominate so that the critical behavior crosses 
over to DP after sufficiently long time. 
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Fig. 46: Simulation of the generalized contact process with n = 2 absorbing states starting from a random 
initial condition. The two different types of inactive domains are shown in black and grey. The active sites 
between the domains are represented by white pixels. 

2. Spontaneous recovery: Active sites turn spontaneously into an inactive state of a 
randomly chosen color. 

3. Boundaries between inactive domains of different colors are free to separate again, 
leaving active sites in between. 

Rules 1 and 2 resemble the usual infection and recovery processes of DP. Rule 3 is new and 
distinguishes different colors. Roughly speaking, this rule ensures that inactive domains 
of different colors cannot stick together irreversibly, rather they will always be separated 
by fluctuating active 'interfaces'. The symmetry under global permutation of the col- 
ors ensures that absorbing domains of different colors compete one another, leading to 
interesting critical behavior. 

Following these descriptive rules, we can introduce a generalized version of the Domany- 
Kinzel cellular automaton (see Sec. |3.2| ) with n absorbing states [263|. It uses the same 
type of lattice and is defined by the conditional probabilities 



(222) 



where k,l = 1, . . . , n and k ^ I. Notice that rule 3 is implemented by transition I^Ii A, 
ensuring that active sites are created between two inactive domains of different colors. 
For n = 1 the above model reduces to the original Domany-Kinzel model. For n = 2 it 
displays a DP2 phase transition at the critical threshold pi^c = P2,c = 0.5673(5). Similarly 
it is possible to define a generalized contact process by the rates 

WAA^Ah = WAA-^hA = l/2n , 
WAik^ikh = wi^A-^hh = 1/2 , ^223) 
yJAIk^AA = Wi^A^AA = A/2 , 







= 1, 


P{A\A,A) = l- 


nP{h\A,A) 


= P2, 


P{A\h,A) 


= P{A\A,h) 


= Pl , 


P{h\h,A) 


= P{h\A,h) 


= 1-Pl 




P{A\hJi) 


= 1, 



wi^i^^i^A = wi^i.^Aii = A/2 , 
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Here the last equation implements rule 3. For n = 1 this model reduces to the usual 
contact process introduced in Sec. For n = 2 the model undergoes a DP2 transition 
at the critical point Ac = 1.592(5). A typical evolution of the generalized contact process 



in 1+1 dimensions is shown in Fig. 4f:. In the active phase A > Ac small inactive islands of 
random color are generated which survive only for a short time. Approaching the phase 
transition their size and lifetime grows while the density of active sites decreases. Notice 
that according to rule 3 a thin film of active sites separates different inactive domains. 
As expected, the numerical estimates of the critical exponents are in agreement with 
the PC exponents ( p21| ). Both models can easily be generalized to higher dimensions. 
However, in higher dimensions the phase transition is presumably characterized by mean 
field behavior. Similarly, increasing the number of absorbing state does not lead to new 
universality classes. Simulations with n > 3 symmetric absorbing states in 1+1 dimensions 
indicate that the system is again described by mean field exponents. 

If the Z2 symmetry of DP2 models is broken by an external field, the critical behavior 



at the transition crosses over to ordinary DP |137]| . Roughly speaking the external field 
favors one of the absorbing states so that pairs of kinks between oppositely oriented inactive 
domains form bound 'dipoles' of a certain size. Interpreting these dipoles as composite 
particles, they recombine and produce a single offspring at certain rates, resembling an 
ordinary DP process on large scales. 

The difference between the PC and DP2 universality classes 

In the DP2 class the 'kinks' between differently colored domains may be interpreted as 
branching-annihilating particles with even number of offspring. Although the number of 
active sites is generally not conserved modulo 2, we may associate with each active island 
between differently colored inactive domains a particle X. Obviously, these particles per- 
form an effective branching-annihilating random walk X — > 3X,2X 0. Therefore, the 
DP2 class and the PC class coincide in 1+1 dimensions. However, it is important to note 
that they are different in higher dimensions. Active sites of PC models in d > 2 dimensions 
can be considered as branching-annihilating walkers, whereas DP2 models describe the dy- 
namics of branching-annihilating interfaces between oppositely oriented inactive domains 
(see Fig. 47). Therefore, the corresponding field theories are expected to be different. A 



field theory for the PC class was presented in Ref. |296, 297|], whereas the development of 



field theories for branching-annihilating interfaces started only recently |24 

In order to understand the difference between PC and DP2, it is helpful to consider two 
other universality classes which also coincide in 1+1 dimensions, namely the annihilation 
process A + A ^ and the voter model [ p52| . The voter model is a two-state model 

with spins Sj = ±1. It evolves by random-sequential updates H >+ + / with equal 

rates. Interpreting H — kinks as particles A, the voter model and the annihilation process 
coincide in 1+1 dimensions. However, in higher dimensions they are different. In fact, 



even their Langevin equations turn out to be different. As shown in Sec. 2.f:, the Langevin 
equation of the annihilation process reads 

dtp{x, t) = -A/)2(x, t) + Z)VV(x, t) + C(x, t) , (224) 
(C(x, t)C(x', t')) = r p2(x, t) - x') 5(t - t') , 

where /9(x, t) represents the coarse-grained density of ^-particles. On the other hand, the 
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Fig. 47: Configurations of an almost critical PC (left) and a DP2 process (right) in d = 2 dimensions. 



Langevin equation of the voter model [253[ is given by 

dtp{x, t) = DV^p{x, t) + C(x, t) (225) 
(C(x, t)C(x', t')) = r /,(x, t) [1 - p(x, t)] d'i^ - x')5(t - t') , 

where p{x, t) £ [0, 1] represents the local orientation of the domain. Obviously, the two 
equations stand for different universality classes. 

DP2 surface critical behavior 

The influence of an absorbing wall in systems with a DP2 transition was studied in 
Ref. |30g[| . Analyzing the generalized Domany-Kinzel model ( |222D with two absorbing 
states in a semi-infinite geometry, it turned out that absorbing and reflective boundary 
conditions play complementary roles. Moreover, since /3 and /3' may be different in the 
DP2 class, one has to introduce two different surface exponents Ps and P'g. For absorbing 
boundary conditions the numerical estimates in a (l+l)-dimensional DP2 process are 

f3s = 1.34(2) , (3i = 2.06(2) . (226) 

For reflecting boundary conditions the two values are simply exchanged. This property 



is related to a duality transformation in parity-conserving processes [310|. It is quite 



remarkable that the two surface exponents seem to obey the scaling relation 

^(/3, + /3^) = i/||-l, {d = l) (227) 
generalizing the conjecture (3^ = i^n — 1 in the case of DP. 



4.3 Activated random walks 

So far we considered nonequilibrium phase transitions where a parameter (e.g. the perco- 
lation probability p) has to be tuned to criticality. Other systems with conserved dynamics 
can be tuned to criticality by varying the particle density in the initial state. An inter- 
esting example is the activated random walk of n particles |23J]. In this model each site 



can be occupied with arbitrarily many particles. Sites with at least two particles are ac- 
tive, i.e., their particles may move independently to randomly selected neighbors. Sites 
with only one particle are frozen. On an infinite lattice, this model has infinitely many 
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absorbing states. The control parameter is the density of particles. For a low density, the 
model quickly evolves into one of the absorbing states, whereas it remains active for high 
particle densities. Near the transition, the stationary density of active sites Pstat scales as 

/0.tat ~ (C - Ccf , C = n/N, (228) 

with the critical point C,c ~ 0.9486 in 1+1 dimensions j23^ ] . However, unlike other models 
with infinitely many absorbing state (see Sec. |3.8| ) the transition does not belong to the 
DP universality class. For example, in 1+1 dimensions the measured exponent /3 = 0.43(1) 
differs significantly from the expected DP value Pdp — 0.277. Obviously, this deviation 
is due to the conservation law. Hence, activated random walks provide a new universality 
class of phase transitions into absorbing states, caused by an additional symmetry, namely 
particle conservation. 



4.4 Absorbing phase transitions and self-organized criticality 

In contrast to ordinary transitions into absorbing states, self-organized critical systems 
exhibit long-range correlations and power laws without being tuned to a certain critical 
point. In some sense the annihilation process A + A ^ (see Sec. ^^ ) can be con- 
sidered as a simple example of self-organized criticality. Without tuning of a parameter 
the annihilation process generates long-range correlations with power-law characteristics. 
Like many other coarsening processes and growth phenomena the driving mechanism is a 
stationary state which is approached but never reached. 

The term 'self-organized criticality' (SOC) refers to a different type of models that 
are attracted to a stationary critical state without being tuned to a critical point. The 



chief examples are the sandpile model and the Bak-Sneppen model |311]. The concept 
of SOC has been used to explain the large variety of power laws observed in nature, for 
example 1// noise, the distribution of earthquakes, and the dynamics of financial markets. 



to name only a few [pl2| . For more than one decade SOC was considered as a quite separate 
field of theoretical statistical physics, being more or less unrelated to conventional phase 
transitions with a tuning parameter. Recently, it was pointed out jSlsHsI^] that SOC 
is in fact closely related to ordinary phase transitions into (infinitely many) absorbing 
states (for a survey see [^38| ) . More precisely, two classes of SOC models have to be 
distinguished. The first class of SOC models, exemplified by the Bak-Sneppen model, 
employs extremal dynamics. In this case the site with an extremal value is selected for 
the next update, i.e., the dynamic rules provide a mechanism of global supervision. In the 
second class, which is represented by the sandpile model, the bulk dynamics is conserved. 
Here a slow driving force competes with the loss of particles at the systems boundaries 
and drives the system to criticality. Hence the process of self-organization is characterized 
by a separation of time scales for avalanches and driving. 

In the latter case, SOC is related to a conventional absorbing-state transition as follows. 



As explained in Ref. | 238 |, any system with conserved local dynamics and a continuous 
absorbing-states transition can be converted into a SOC model by (1) adding a process for 
increasing the density in infinitesimal steps, and (2) implementing a process for decreasing 
the density infinitesimally while the system is active. For example, let us consider the 
activated random walk. Adding a process for random deposition of particles and a process 
for loss of particles during avalanches at the boundaries, we obtain the so-called Manna 
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sandpile model |316] which is known to exhibit SOC. Using a deterministic variant of the 
same model, one obtains the famous Bak-Tang-Wiesenfeld model [^]. However, in order 
to create a SOC counterpart for directed percolation, a slightly different recipe has to be 



used, as demonstrated in Ref. [314, 317 1. 



4.5 The annilation/fission process 



As shown in Sec. 2.6, the influence of fluctuations in (l+l)-dimensional systems can be 



very different. In the annihilation process, for example, the particles become anticorre- 



lated, i.e., they try to be far away from each other. As shown in Ref. [43|, this type of 
fluctuations corresponds to 'imaginary' noise in the Langevin equation. In a DP process, 
however, the particles are highly correlated and the noise turns out to be real. Three 
years ago Howard and Tauber posed the question whether it is possible to interpolate 
between real and imaginary noise. As a prototype of such a transition, they considered 
the annihilation/fission process 

2A^3A, 2A^ (229) 

with diffusion of single particles. Note that this is a binary process, i.e., at least two par- 
ticles are required to meet at the same (or at neighboring) sites in order to self-destruct 
or to create offspring. Moreover, there are no exceptional symmetries such as parity con- 
servation on the microscopic level. The model exhibits a nonequilibrium phase transition 
between an active phase and two non-symmetric absorbing states, namely the empty lat- 
tice and the state with only one diffusing particle. Performing a field-theoretic Howard 
and Tauber argued that this transition should belong to an independent yet unknown 
universality class of phase transitions. 

The annihilation/fission process may be interpreted as a pair contact process plus 
diffusion of single particles. As a consequence, the static background shown in Fig. 
begins to fluctuate. Since the decay of the particle density in the subcritical phase is 
governed by the annihilation process 2A — > 0, it is natural to expect that the transition 
does not belong to the DP class. Similarly, the possibility of DP2 critical behavior seems 
to be unlikely since there is no parity conservation or Z2-symmetry in the system [ 163| , pl8 |. 



Thus, the annihilation/fission process may represent a new type of nonequilibrium critical 
phenomenon which has not yet been studied before. 

4.6 The Lipowski model 

A particularly interesting model with a non-DP transition in d = 2 spatial dimensions has 



been introduced recently by Lipowski |31£]. It has infinitely many absorbing states and is 



defined by extremely simple dynamic rules. 

The model is defined as follows. The state of the system is specified by real numbers 
Uij G (0,1) which reside on the bonds of a d-dimensional square lattice (see Fig. 48). 



Initially all these variables are randomly distributed between and 1. The model evolves 
by random sequential updates depending on a control parameter p which plays the role of 
a percolation probability. For each update attempt a site i is selected at random. Its local 
'field' hi is defined by the sum of the y-variables of the four connecting bonds. li hi > p 
the four variables are replaced by random numbers drawn from a fiat distribution between 
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Fig. 48: The Lipowski model. Each bond of a square lattice carries a weight between zero and one. A 
plaquette of four bonds (hollow circles) is randomly selected. If the sum of the four weights exceeds a 
certain threshold, they are replaced by four independent random numbers between zero and one. 

and 1. Otherwise the update attempt is abandoned. As usual in models with random- 
sequential dynamics, each update attempt corresponds to a time increment of At = 
where A'' is the total number of sites. A site is called 'active' if it is susceptible for a 
replacement, i.e. hi > p. Measuring the density of active sites in a numerical simulation, 
the model displays a continuous phase transition between a fluctuating active and a frozen 
phase. Clearly, the model has infinitely many absorbing states. 

The Lipowski model poses a puzzle: In one dimension the static exponent P coincides 
with the DP exponent /? ~ 0.276. In two dimensions, however, the static exponent P is 
found to be much smaller than the expected value 0.58. Surprisingly, the measured value 
seems to coincide with the DP value in one dimension. Lipowski argued that the model 
should provide a mechanism for dimensional reduction, placing the two-dimensional critical 
phenomenon into a one-dimensional universality class. However, the observed coincidence 
may well be accidental. Moreover, it is not obvious how such a mechanism should work. It 
is therefore an interesting open question whether dimensional reduction can be observed 
in stochastic lattice models. 
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5 Damage spreading 



One of the central problems of dynamic system theory is the dependence of the system's 
temporal evolution on the initial conditions. It is well known that nonlinear systems with 
deterministic dynamics may be extremely sensitive to small perturbations of the initial 
state. Even in simple systems such as in a periodically driven pendulum, a small variation 
of the initial parameters can change the entire trajectory completely. If the distance 
between two infinitesimally different trajectories diverges during the temporal evolution, 
a dynamic system is said to exhibit chaotic behavior. 

The notion of chaos has been introduced in the context of deterministic systems where 
a trajectory is uniquely determined by the initial condition. In random processes, however, 
trajectories are not uniquely determined; the time evolution of such a system is not repro- 
ducible and therefore the usual definition of chaotic behavior does not apply. Nevertheless, 
one may pose the question how the system responds to changes in the initial condition. 
It would be interesting to know, for example, how the biological evolution on earth would 
have been affected if the initial conditions were slightly different. In order to address 



this question, Kauffman introduced the concept of damage spreading (DS) [|320 , 321 1. In 



damage spreading simulations 1 322 ] two copies (replicas) S, S' of a stochastic model evolve 
simultaneously. Initially the two copies differ only on a small number of sites. This differ- 
ence is considered as a small perturbation (damage) in one of the two systems. Moreover, 
it is assumed that the replicas evolve under identical realizations of thermal noise, i.e., 
both copies use the same sequence of random numbers in the simulation. If the number 
of sites in different states, the so-called Hamming distance 

A(t) = Mt) = E 1 - Ut),s[it) (230) 

i i 

does not go to zero in the long-time limit, damage is said to spread, indicating high 
sensitivity with respect to the initial condition. Otherwise, if A{t) vanishes, damage is 
said to heal, indicating a weak influence of the initial condition. In order to get statistically 
meaningful results, the Hamming distance has to be averaged over many realizations of 
randomness. 



Damage spreading first appeared in physics literature in the mid eighties 1 322 -327] 
and attracted considerable interest and attention. The main reason behind this initial 
enthusiasm was the hope that damage may spread in some regions of a system's parame- 
ter space and disappear elsewhere, indicating the existence of chaotic and regular phases 
in stochastic systems. The initial enthusiasm abated during subsequent years, the main 
reason being an apparent lack of an objective measure whether damage does or does not 
spread in a given system. More precisely, it was realized that the location of the phase 
boundaries may depend on details of the algorithmic implementation. However, if spread- 
ing or healing of damage indicated some intrinsic property of the system, one would not 
expect the result to depend on details of the algorithm used to generate its dynamics. 
Meanwhile DS has been applied to a large variety of models (see Table In view of 
the vast literature on DS simulations, it is therefore necessary to carefully analyze the 
conceptual problems of this technique. 

In the following we discuss several aspects of DS. First we present a simple example in 
order to explain how a DS simulation depends on the algorithmic implementation. From a 
more mathematical point of view, this phenomenon can also be understood by analyzing 
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Equilibrium models: 


Ref. 


Ising ferromagnet 
Heisenberg model 
XY model 

Potts and Ashkin-Teller models 
discrete A^- vector ferromagnets 
spin glasses 

Ising model with microcanonical constraints 

diluted Ising model 

Ising model in quenched random fields 

Frustrated Ising systems 

Layered Ising systems 


1 
1 
1 
1 

1 
1 

1 
1 

1 
1 
1 


525 
532 
53S 
534 
537 

m 

546 
547 
54£ 
55C 
551 


,326 

-336 
-345 
,348 


,|328-331| 


Nonequilibrium models: 




Kauffman model 
The game of life 
ZGB model 

Ohta-Jasnov-Kawasaki and Ginzburg-Landau model 

Irreversible reaction-diffusion processes 

Models for surface reactions 

Restricted solid on solid growth models 

Models of the PC class 

SOC models 

Traveling salesman problem 
Boolean random networks 


1 
1 
1 
1 
1 

1 
1 
1 

1 
1 
1 


52C 
554 
556 
557 

56C 
561 
562 
56S 
564 
565 


321 
,355 

359 


,|352,^5S1 



Table 3: Some applications of damage spreading. 

the joint master equation of the two replicas. Because of their algorithmic dependence, 
DS simulations are ambiguous and cannot be used as a criterion for 'chaotic' and regular 
phases. However, to some extent the ambiguity of DS can be overcome by an algorithm- 
independent definition of DS phases. In this approach the entire family of physically 
legitimate algorithms for a given dynamic system is considered as a whole. Furthermore, 
we summarize what is known about the universal properties of DS transitions. Finally we 
discuss several applications of DS simulations. 



5.1 Damage spreading phases 



Algorithmic dependence of damage spreading simulations 

The conceptual problem of DS was first discovered in the Domany-Kinzel cellular automa- 
ton (cf. Fig. ^). Martins et al. [366| observed that in a certain region of the active phase 
damage spreads and heals elsewhere. Subsequently several other authors determined the 
boundary of this region with increasing accuracy |114, 367-p6£ 



Independently, mean- 
field type approximations of varying complexity confirmed the existence of this 'chaotic 
phase' [ p68 , [370 - 372 1 . Its boundary, however, was shown to depend on the manner in which 
the dynamic procedure of the DK model is implemented on a computer |370 , 37l|] , while the 
evolution of a single replica is completely insensitive to the algorithmic implementation. 



This prompted Grassberger 369 1 to observe that 



"it is misleading to speak of different phases in the DK automaton, ...instead 
these are different phases for very specific algorithms for simulating pairs of 
such automata" . 
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equilibrium 
model 

Ising model 
Potts model 
Ashkin-Teller 
N-vector 




dynamic 

rule 

f ^ 
heat bath 

Metropolis 

Kawasaki 

Q2R 

Wolf 




algorithmic 
implementation 

spin-orienting (heat bath) 
spin-flipping (Glauber) 
mixtures of Glauber / heat bath 



Fig. 49: The thermodynamic ensemble of an equilibrium model (left) may be generated by different 
dynamic procedures (center) that are characterized by specific transition probabilities. Each of these 
dynamic procedures in turn can be realized by difi'erent algorithms (right). These algorithms are fully 
equivalent and cannot be distinguished by an observer of a single system. However, their damage spreading 
properties turn out to be different (see text). 



To understand the algorithmic dependence of DS simulations, let us consider a much 
simpler system, namely the Ising model. In this context it is important to note that 



there are different levels of variety in stochastic lattice models (see Fig. 4£ ) . As explained 
in Sec. p.lO] , the equilibrium ensemble of the Ising model can be generated by various 
different dynamic rules. For example, heat bath and Metropolis dynamics represent two 
different dynamic systems which have the same stationary state. Initially it was hoped 
that DS would not depend on the intrinsic dynamics, allowing regular and chaotic phases 
to be identified as equilibrium properties [ |325| -327|. However, later it was realized that 
different dynamic procedures (such as Glauber versus Metropolis 1 32S| , 331 1 , Q2R [l329|l , or 



Kawasaki dynamics [330|) exhibit different DS properties. Yet, this is not surprising since 
the dynamic rules, although generating the same equilibrium ensemble, represent different 
dynamic systems. Since DS is a dynamic phenomenon it is quite natural that it depends 
on the dynamics under consideration. Similarly, it is not surprising that DS depends on 
the type of updates [373, p74 | and the interaction range 375 1. 

The real conceptual problem of DS occurs at the second level in Fig. On a com- 
puter each dynamic rule can be realized by several algorithms. Regarding a single system, 
these algorithms are fully equivalent and cannot be distinguished. However, in DS simula- 
tions they lead to different results. To understand this apparent paradox, let us consider 
the Ising model with spin-orienting (standard heat bath) and spin-flipping (Glauber) dy- 
namics. Both algorithmic prescriptions represent the same dynamic rule which mimics 
the contact of the Ising model with a thermal reservoir by means of local spin dynamics. 
More precisely, an observer of a single system who analyzes the trajectories would be 
unable to decide whether its temporal evolution was generated by standard heat bath or 
Glauber dynamics. This can already be seen in the example of a single-spin Ising model 
at infinite temperature. Since in this case the spin (T(i) = ±1 changes randomly in time, 
the transition rates are simply given by W-\^j^\ = = 1. These transition proba- 

bilities define our dynamic system. However, this system can be realized by two different 
algorithms, namely by spin-orienting updates (standard heat bath dynamics) 

z=rnd(0,l) ; 

if (z<0.5) a(t^dt) = \\ else a(t^ dt)=-\; 
and by spin-flipping updates (Glauber dynamics) 
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z=rnd(0,l) ; 

if (z<0.5) a{t + dt)=-a{t); else a {t + dt)=a{t) . 

It is obvious that both procedures are fully equivalent on a single replica, i.e., an observer 
would be unable to decide which of the two algorithms has been used. The difference 
between the two algorithms may become evident only if we observe the evolution of two 
replicas S and S' of the system in a DS simulation: For spin-orienting updates an initial 
'damage' a{0) = — (j(0)' heals immediately while it is preserved when spin-flipping updates 
are used. 



A similar algorithmic dependence can be observed in the full Ising model (| 
temperature. Defining the transition probability (cf. Eq. (p9[)) 



^hi(t) _|_ Q-hi{t) 

we may express the update rules of heat bath dynamics by 

ai(t + l) = sign[pi{t) - Zi{t)] , 



at finite 



(231) 



(232) 



where Zi{t) G (0,1) are equally distributed random numbers. The corresponding update 
rule for Glauber dynamics is given by 



fJi{t + l) 



+sign[pi{t) - z] 
-sign[l -pi{t) 



if ai{t) 



+1 
-1 



(233) 



It is easy to verify that for given {crj_i(t), crj(t), (Tj+i(t)}, the probability to get crj(t + 1) = 
-|-1 is the same in both cases. Hence, by observing the temporal evolution of a single Ising 
system, one cannot tell which of the two methods was used to generate the trajectory in 
configuration space. The two algorithms can only be distinguished when two copies are 
simulated in parallel, i.e., by studying damage spreading. 

Investigating the two-dimensional Ising model with Glauber dynamics, Stanley et 



al. 1 326] and Mariz et al. [328| found that damage spreads in the disordered phase T > Tc 
and heals elsewhere. Performing more precise simulations, Grassberger p76| ] realized that 
the DS transition occurs slightly below Tc- This observation was also be supported by a 



mean field theory [377|. Moreover, the critical point of the DS transition was found to vary 



continuously when mixtures of Glauber and heat bath dynamics are used [375]. Very simi- 
lar properties are observed in the three-dimensional Glauber model [p76 ,379 ~ p81 |. Turning 
to the Ising model with spin orienting (standard heat bath) dynamics, it is possible to 



prove that damage does not spread at any temperature in any dimension [327, 328 [. 



The master equation of damage spreading simulations 

In order to understand the algorithmic dependence of DS simulations from a more funda- 
mental point of view, let us consider two copies Si and 52 of an arbitrary nonequilibrium 
system with asynchronous dynamics (random-sequential updates). Each of the two sys- 
tems evolves according to a master equation with a given Liouville operator C. The total 
system S = (51,52) constitutes a new nonequilibrium system that evolves according to 
a joint master equation with a certain Liouville operator ^A. If both systems were using 
independent random numbers, M would be given by the tensor product C. However, 
in DS simulations the use of the same sequence of random numbers leads to nontrivial 
correlations between 5i and 52. 
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According to the definition of damage spreading, the Liouvihe operator M. of the total 
system (5*1, ^2) is restricted by certain physical constraints. On the one hand, each replica 
is required to evolve according to its own natural dynamics. Hence, by integrating out 
the degrees of freedom of one of the replicas we obtain the Liouville operator of the other 
system: 

J2{si,S2\M\s[,s'2) = {s2\C\s'2) ioT all S2, s[, s'2 , (234) 

si 

''^^{si, S2\M.\s[, s'2) = {si\C\s[) for all Si, s'l, S2 • (235) 



S2 



These restrictions already imply probability conservation for the total system. On the 
other hand, the trajectories of the two replicas, once they have reached the same state (no 
damage), have to be identical: 

(si,S2|A^|s',s') = (si|£|sVsi,s2 . (236) 

The ambiguity of damage spreading simulations is due to the fact that the restrictions 



( 234 )-( p36D do not fully determine the Liouville operator A4 of the total system. This can 
easily be verified by counting the degrees of freedom. For a system with n configurations 
the restrictions determine less than 3n^ of the matrix elements of A4 so that damage 
spreading is ambiguous for n > 3. But even in the case of a single-spin Ising model with 
n = 2 states one can show that only 14 of 24 equations are linearly independent so that two 
out of 16 matrix elements of A4 can be chosen freely. The remaining degrees of freedom are 
the origin of the algorithmic dependence of DS simulations. A similar ambiguity occurs 
in the joint transfer matrix of models with parallel dynamics ||382(|. 



An algorithm-independent definition of damage spreading 

In order to overcome these conceptual difficulties, let us consider the entire family of 
dynamic procedures consistent with certain physically dictated constraints |382|. Then 
for any particular system one of the three possibilities may hold: 



1. Damage spreads for any member of the family of dynamic procedures. 

2. Damage heals for any member of this family. 

3. Damage spreads for a subset of the possible dynamic procedures, and heals for the 
complementing subset. 



Obviously this allows us to classify damage spreading properties in an algorithm-indepen- 
dent manner. To this end we must, however, consider simultaneously the entire set of 
possible algorithms (dynamic procedures) which are consistent with the physics of the 
model under consideration. This set of physically 'legitimate' algorithms may be defined 
by certain restrictions that are dictated by the dynamics of the single evolving system: 
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Fig. 50: Algorithm-independent damage spreading phases in the Domany-Kinzel model. 1: Damage never 
spreads. 2: Damage may spread, depending on the algorithm. 3: Damage always spreads. 



a) Definition of DS: The dynamic rules for the evolution of the pair of replicas are such 
that a single replica evolves according to its 'natural' dynamics. Once both replicas 
have reached the same configuration, their temporal evolution will be identical. 

b) Interaction range: The transition probabilities for the combined system at site i may 
depend only on those sites that affect the evolution of site i under the dynamic rules 
of a single system. 

c) Symmetry: The rules that govern the evolution for the pair of systems do not break 
any of the symmetries of the single-replica dynamics. 



The first restriction is simply a verbal formulation of Eqs. ( p34| )-( p36| ). The second con- 
dition tells us that the interaction range in the combined system of two replicas must 
not exceed the interaction range of a single system. The third rule implies, for example, 
that if there is a left-right symmetry in the evolution of a single system, the same must 
hold for the pair of replicas. It can be shown that these conditions suffice to unambigu- 
ously determine a set of physically 'legitimate' algorithms. This was demonstrated for the 



one-dimensional Domany-Kinzel model [382] for which the phase diagram in Fig. 5C was 
found. 

Clearly, the subjectivity in defining DS phases, as described before, has now been 
shifted to selecting the restrictions defining which DS procedure is 'legitimate'. However, 
the specification of such a family by means of physically motivated criteria appears to be 
less arbitrary than choosing, at random, one out of a continuum of physically equivalent 
update procedures. It should also be emphasized that the algorithm-independent definition 
of DS phases does not mean that DS is reflected in the dynamic behavior of a single system, 
so that Grassberger's observation still holds: DS is a property of a pair of stochastic 
systems. 



106 



replica 1 



replica 2 



damage 



Fig. 51: Damage spreading in the (l+l)-dimensional Domany-Kinzel cellular automaton evolving in the 
stationary active state at pi = 0, P2 = 0.313. An initial damage is introduced at the lattice site in the 
center. 



5.2 Universality of damage spreading transitions 
The DP conjecture for damage spreading 



As can be seen in Fig. |51|, spreading of damage is in many respects similar to spreading 
of activity in a DP process. As in DP damage spreads to nearest neighbors and heals 
spontaneously. Once both copies have reached the same configuration (no damage), their 
evolution will be identical, i.e., they are confined to some 'absorbing' subspace from where 
they cannot escape. In contrast to DP the spreading process depends crucially on the 
actual evolution of the two replicas, providing a fluctuating background in which the 
spreading process takes place. Nevertheless DS follows the same spirit as the spreading of 
activity in a DP model. This observation led Grassberger to the conjecture that damage 
spreading transitions belong generically to the directed percolation universality class p69|] . 



So far analytical support for this conjecture came from approximate mean-field argu- 
ments [371] and an exact statement by Kohring and Schreckenberg p70[|, who noted that 



on the P2 = line the dynamics of damage spreading in the DK automaton is precisely 
identical to the evolution of the DK automaton itself, hence on this line DS is trivially in 
the DP universality class. To prove this statement, consider two replicas of S and S' of a 
DK automaton evolving according to Eq. ( |83| ) with equal random numbers Zi[t) = -z^(t). 
The probabilities Pni^i = s[_i,s[_^^) to generate damage at site i are listed 

in Table ^. For P2 = they may be expressed as 

= Si+i and s^_i / s^+i , 
Pd{\ = l|si-i,Si+i; s-_i,s-+i) = { pi if Si-i / Sj+i and s'^_-^ = s'^^-^ , (237) 




or equivalently 

Pd{\ = l|si_i, Sj+i; S'+i) = pi(l - (5a,_i,a,+i) • (238) 

Therefore, the damage variables Aj(t) evolves precisely according to the probabilistic 
rules of a single DK automaton. Hence for P2 = the DS transition belongs to the DP 
universality class. This mapping of DS to DP was later extended to other regions in the 



phase diagram of the DK model [382|. Although such an exact mapping is usually not 



available for other models, various numerical simulations show that most DS transitions 
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Table 4: Probabilities PoiAi — Si+i; s^_i, s^+i) for the generation of damage in the DK model. 



are indeed characterized by DP exponents, supporting Grassberger's conjecture. The 
same apphes to deterministic cellular automata with chaotic behavior if a small noise is 



added p3 | 



Different DS properties may be expected in models with cluster dynamics [384[ . How- 



ever, as pointed out in Refs. [385, 386| ], it is difficult to extend the definition of 'using 



the same random numbers' to cluster algorithms since the number of clusters may differ 
on both replicas. This difficulty can be overcome by introducing a random background 
field [ |387]] . Assigning a local random number to each lattice site the new orientation of a 
cluster depends on whether the sum of its random numbers positive or negative. Although 
cluster algorithms involve long-range correlations, DS transitions still belong to the DP 
universality class unless they coincide with the thermodynamic phase transition of the 
Ising system. 

DS transitions vi^ith non-DP behavior 

The critical properties of a DS transition are expected to change if one of the four con- 
ditions of the DP conjecture is violated. For example, DS transitions in models with 
frozen randomness (such as the Kauffman model [ ^20| , |32l| ) do not belong to the DP class. 
Non-DP behavior is also expected when the DS order parameter exhibits an additional 
Z2 symmetry (cf. Sec. [1.2D . In the context of damage spreading it is important to note 
that such a symmetry should be a property of the DS order parameter, i.e., the Hamming 
distance. For example, the Z2 symmetry of Ising systems is not sufficient - inverting all 
spins in both replicas does not change the Hamming distance between the two configura- 
tions. Similarly, models with a non-DP transition do not automatically exhibit non-DP 
damage spreading. For example, Grassbergers cellular automaton A, which has a DP2 
phase transition, displays an ordinary DS transition belonging to DP [p62|| . 

Non-DP behavior at the DS transition can be observed in systems with a symmetry 
between two 'absorbing states' of Aj(t), one without damage A = and the other with 
full damage A = 1. The simplest example of such a symmetry is the Ising model with 
Glauber dynamics. In fact, if both replicas are in opposite states (full damage), they will 
always evolve through opposite configurations. Unfortunately, there is no DS transition 
in the one-dimensional Glauber model. However, by exploiting the algorithmic freedom, 
it is possible to construct a modified dynamic rule which exhibits a DP transition |l388|| . 
This rule depends on a parameter A and is defined by 

+sign(pi(t) -z) ify = l . . 

^^^' + ^^-1 -sign(l-p.(t)-z) ify = -l ' ™ 



1 + (Ji_i(t)cJi+i(t) ) + ( 1 - (Ti_i(t)c7i+i(t) )sign(A - z 
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where 2, z G (0, 1) are two independent random numbers. On a single replica this update 
rule is fully equivalent to Glauber and heat bath dynamics for all < A < 1. However, the 
effective rate for spreading of damage depends on A. For fixed temperature J/ksT = 0.25 
a DS transition with non-DP exponents occurs at the critical value A* = 0.82(1). This 
example demonstrates that additional symmetries of the DS order parameter may lead to 
non-DP behavior at the DS transition. A similar situation emerges in the nonequilibrium 
Ising model introduced by Menyhard [p62| . 



5.3 Applications of damage spreading 

Measurement of critical exponents in equilibrium models 

Damage spreading simulations can also be used to determine certain static and dynamic 
properties of systems at thermal equilibrium. This application of DS was first demon- 



strated by Coniglio et al. [389| who showed that there exists an exact relation between the 
Hamming distance A and certain correlation functions. The essential idea is to consider 
two copies S, S' of an Ising model with spins cr^ = ±1 and introducing a small damage 
by keeping a single spin in one of the systems fixed during the whole temporal evolution, 
say (Tq = —1. Since this perturbation breaks the symmetry between the two copies, it is 
important to distinguish two different types of damage at site i, namely cjj = l,a[ = — 1 
and fjj = —l,cr'- = 1. The probabilities of finding these different types of damage in the 
stationary state is given by 

d+- = ((1 - s^)s'^ , d-+ = {si{l - Si)') , (240) 

where Si = \{cFi — 1) (see also Ref. i322[ ) . Notice that these quantities can be understood 
as two-point correlation functions between the two copies. However, their difference 

T, = dt- -dT+ = {s,)-{s',) (241) 

is a combination of one-point functions, i.e., by taking a certain combination of the dam- 
age probabilities one obtains quantities that describe the properties of a single system. 
Therefore, such quantities do not depend on the algorithmic implementation. In fact, 
using detailed balance and ergodicity one can prove that 

^ - (242) 



' 2(1 -m) ' 

where Cq^j is the two-point correlation function and m the magnetization of the Ising 
model at thermal equilibrium. This relation is exact and does not depend on the specific 
algorithmic implementation used in the simulation. Moreover, it can be shown by mono- 
tonicity arguments that for standard heat bath dynamics the total damage A is related 
to the static susceptibility x by 

X = 2{l - m)Y,^i . (243) 

i 

Thus, damage spreading simulations can be used to determine the static exponents (3 and 



and the critical temperature of the Ising model (see Refs. [39C- ^9^ ). Similar relations 
between Hamming distance and correlation functions were found in certain lattice models 
with absorbing states |l393|| . 
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Fig. 52: Identification of domain walls. Left: Tlie two-dimensional Ising model evolving by h eat bath 
dynamics. Right: C orres ponding domain walls detected by the observable Ai defined in Eq. (244). Figure 
reprinted from Ref. 



DS is also used as a tool for accurate measurements of the dynamic exponent z at the 
phase transition of equilibrium systems, provided that the thermodynamic transition and 
the DS transition coincide [|394 , 395(| . In this case the spreading exponent z is not given by 
the DP exponent vs^jv^^ instead it is expected to coincide with the dynamic exponent z of 
the model under consideration. The knowledge of z is important in order to estimate the 
critical slowing down of given dynamic system. In a series of papers the dynamic exponent 
of the Ising model with heat bath dynamics has been measured in two ||396| - f^00|| and three 
dimensions [397, 39^, 401, 402 1. After an initial controversy, Grassberger |403| compared 
several simulation methods and found the estimates z = 2.172(6) in two and z = 2.032(4) 
in three dimensions. Later refined simulations confirmed these results |404|. 



Identification of domain walls in coarsening systems 

Damage spreading techniques can also be used to identify domain walls in coarsening 



systems. Coarsening phenomena [ 405 | occur in various dynamic systems as, for example, 
in the Ising model with Glauber dynamics. In the ordering phase, starting with random 
initial conditions, patterns of ferromagnetic domains are formed whose typical size grows 
with time as t^/"^. For zero temperature, these domains are completely ordered and the 
domain walls can be identified as bonds between oppositely oriented spins. For nonzero 
temperature, however, it is difficult to define domain walls as one has to distinguish be- 
tween 'true' domains and islands of the minority phase generated by thermal fluctuations. 
In order to identify coarsening domains for < T < Tc, Derrida [|217| proposed to compare 
two identical copies evolving under the same thermal noise. One copy starts with random 
initial conditions and begins to coarsen, whereas the other copy starts from a fully magne- 
tized state and remains ordered as time evolves. It is assumed that all spin flips occurring 
in ordered replica can be regarded as thermal fluctuations. Therefore, when a spin flip 
occurs simultaneously in both replicas, it can be considered as a thermal fluctuation, while 
it is a signature the coarsening process otherwise. This method was used to determine the 



fraction of persistent spins of the Ising model |207] as a function of time, confirming that 
the persistence exponent does not change for < T < Tc. 

Since Derrida's method allows only one type of domains to be detected, it is impossible 
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to identify the precise location of domain walls. This can be overcome by comparing three 



replicas instead of two |219|. As before, the first copy starts with random initial conditions 
and serves as the master copy in which the coarsening process takes place. The other 
copies start from fully ordered initial conditions with positive and negative magnetization, 
respectively. Fluctuations in the first copy, which do not occur in the other two replicas, 
indicate the presence of a domain wall. More precisely, domain walls may be detected by 
the observable 

where the superscripts denote the replica. As shown in Fig. |5^, this method works sur- 
prisingly well. However, it turns out that the appearance depends on the algorithmic 
implementation, i.e., Glauber and heat bath dynamics leads to different results. This is 
not surprising as the method relies on the same ideas as damage spreading. 



5.4 Damage spreading and experiments 

In this Section we have seen that the concept of damage spreading depends on the algorith- 
mic implementation and therefore lacks a well-defined meaning. The suggested algorithm- 
independent definition of DS resolves these difficulties only partly, since it is based on 
certain (physically motivated) ad hoc assumptions. Therefore, DS does not provide a 
strict definition of 'chaotic' and 'regular' phases in stochastic systems. Nevertheless DS 
has been useful to estimate critical exponents of certain nonequilibrium phase transitions 
and to identify domain walls in coarsening system at nonzero temperature. Concerning 
the question whether DS is relevant for experiments, there is a clear answer: DS is an 
artificial concept which does not exist in nature. In particular, there is no meaning of 
'using the same sequence of random numbers' in an experimental system. DS is rather a 
simulation technique for a pair of systems, taking advantage of our ability to use determin- 
istic pseudo random number generators. It is indeed not surprising that such a concept, 
for which there is no correspondence in nature, expresses its incompleteness by certain 
inconsistencies. 



Ill 



6 Interface growth 



A further important field of statistical physics is the study of crystal growth and transitions 
between different morphologies of moving interfaces |406| . During the last two decades 
there has been an enormous progress in the understanding of growth processes (for reviews 
see e.g. Refs. |33, 407-409|). In this Section we will focus on certain classes of depinning 
transitions which are closely related to nonequilibrium phase transitions into absorbing 
states. 



6.1 Roughening transitions — a brief introduction 

Models of growing interfaces may be realized either on a discrete lattice or by continuum 
equations. Discrete solid-on-solid (SOS) models are usually defined on a d-dimensional 
square lattice of N = L'^ sites associated with integer variables hi marking the actual 
height of the interface at site i. Clearly, this description does not allow for overhangs 
of the interface. The set of all heights {hi} determines the state of the system which 
evolves according to certain stochastic rules for adsorption and desorption. In restricted 
solid-on-solid models (RSOS) the dynamic rules satisfy the additional constraint 

\h,-h,+i\<l, (245) 

i.e., adjacent sites may differ by at most one height step. To characterize the evolution of 
the interface, it is useful to introduce the mean height h{t) and the width w{t): 

i i i 

A surface is called 'smooth' if the heights hi are correlated over arbitrarily large distances. 
Otherwise, if distant heights become uncorrelated, the surface is said to be rough. A rough 
surface is typically characterized by a diverging width when the system size is taken to 
infinity. 

In many cases growing interfaces exhibit simple scaling laws. In a finite system, starting 
from a flat interface, the width first increases algebraically as w ^ t"' , where 7 > is the 
growth exponent of the systemf^. In the long-time limit, when the correlation length 
reaches the system size, the width saturates at some constant value Wgat ~ L", where 
a is the roughness exponent. This type of scaling behavior is known as Family-Vicsek 
scaling | |410| ], as described by the scaling form 

wiN,t)^N^fit/N') , (247) 

where z = a/7 is the dynamic exponent of the model. f{u) is a universal scaling function 
with the asymptotic behavior f{u) ~ for n — > and f{u) = const for u ^ 00. Notice 
that this scaling form is invariant under rescaling 

X ^ Ax , t^AH, /i(x, t) A"/i(Ax, AH) , (248) 



^^In most textbooks the growth exponent is denoted by 13. Here we use a different symbol in order to 
distinguish this exponent from the density exponent of DP. 
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where A is a dilatation parameter. Notice that in contrast to the order parameter of a 
spreading process (^), the fluctuations in the heights increase under rescaling. 

The critical exponents 7 and z are used to categorize various universality classes of 
roughening interfaces. Typically each of these universality classes is characterized by 
certain symmetry properties of the system and may be associated with a specific stochastic 
differential equation. This equation describes the growth of a continuous height field /i(x, t) 
and consists of the most relevant operators under rescaling ( |248| ) that are consistent with 
the symmetries of the system. For example, postulating the symmetries 



1. translational invariance in space x — > x + Ax, 

2. translational invariance in time t — > t + At, 

3. translational invariance in height direction h ^ h + Ah, 

4. reflection invariance in space x — x, 

5. up /down symmetry h ^ —h, 



one is led to the Edwards-Wilkinson (EW) universality class |411], as described by the 
equation 

= V + aV'h{^, t) + C(x, t) , (249) 

where v denotes the mean velocity and a the surface tension. C(x, t) is a zero-average 
Gaussian noise field with variance 

(C(x, t)C(x', t')) = 2D6'^{x - x')5{t - t') (250) 

taking the stochastic nature of deposition into account. This equation is linear and thus 
exactly solvable. Scaling invariance ( |24§| ) implies that the critical exponents are given by 

a = 1 - d/2 , 7 = 1/2 - (i/4 , z = 2 . (EW) (251) 

Edwards- Wilkinson growth processes are invariant under up/down reflection of the in- 
terface h — > —h. However, if atoms are adsorbed from a gas phase above the interface 
there is no particular reason for the system to be up/down symmetric. In that case the 
above equation has to be extended by the most relevant term that breaks the up/down 
symmetry, leading to the Kardar-Parisi-Zhang (KPZ) equation |18f:, 412| 



^M^ill = y + ^v2/i(x, t) + A(V/i(x, t)f + C(x, t) . (252) 
ot 

For a one-dimensional interface the critical exponents of the KPZ universality class are 
given by 

a = 1/2, 7=1/3, z = 3/2. (KPZ in Id) (253) 

whereas in d >2 dimensions only numerical estimates are known (see Ref. [^). 

It is particularly interesting to study roughening transitions between a smooth and 
a rough phase. A roughening transition is usually accompanied by a diverging spatial 
correlation length In the smooth phase this correlation length provides a typical scale 
below which the interface appears to be rough. On larger scales, however, the interface 
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turns out to be smooth. Approaching the roughening transition the correlation length 
diverges whereby the entire interface becomes rough. One of the simplest models displaying 
a roughening transition is the two-dimensional discrete Gaussian SOS model [^. Another 
important example is the KPZ equation which exhibits a roughening transition in d > 2 
spatial dimensions. 

In the following we will focus on certain growth models with depinning transitions. 
In particular we will discuss depinning transitions in random media, polynuclear growth 
processes, and solid-on-solid growth processes with evaporation at the edges of plateaus. In 
the pinned phase of these models the interface is smooth and does not propagate. Varying 
a control parameter the interface undergoes a depinning transition; it starts moving and 
evolves into a rough state. As we will see below, various depinning transitions are closely 
related to phase transitions into absorbing states. 



6.2 Depinning transitions of driven interfaces 



An interesting class of depinning transitions can be observed in experiments of driven 
interfaces in random media |3^. In these experiments a liquid is pumped through a 
porous medium. If the driving force F is sufficiently low the liquid cannot move through 
the medium since the air/liquid interface is pinned at certain pores. Above a critical 
threshold, however, the interface starts moving through the medium with an average 
velocity v. Close to the transition, v is expected to scale as 

vr^iF- F,f , (254) 

where 6 is the velocity exponent. One of the first experiments in 1-1-1 dimensions was 
performed by Buldyrev et al., who studied the wetting of paper in a basin filled with 



suspensions of ink or coffee |412]. Here the driving force F is a result of capillary forces 
competing with the total weight of the absorbed suspension. Consequently, the interface 
becomes pinned at a certain height where F ~ Fc. Once the interface has stopped, the 
width should scale as 

w{i,t)^rfit/f), (255) 

where £ is a box size and a the roughening exponent. Measuring the interface width 
Buldyrev et al. found the roughness exponent a = 0.63(4). In various other experiments 
the values are scattered between 0.6 and 1.25. This is surprising since the Kardar-Parisi- 
Zhang (KPZ) class |186| predicts the much smaller value a = 1/2. 



It is believed that the large values of a are due to inhomogeneities of the porous 
medium. Because of these inhomogeneities, the interface does not propagate uniformly by 
local fluctuations as in the KPZ equation, it rather propagates by avalanches. In the lit- 
erature two universality classes for this type of interfacial growth have been proposed. In 
case of Hnear growth the interface should be described by a random field Ising model ||414|| , 
leading to the exponents a = 1, -2 = 4/3, and 9 = 1/3 in 1-1-1 dimensions. In the presence 
of a KPZ-type nonlinearity, however, the roughening process should exhibit a depinning 
transition which is related to DP [ 415| , 416|. The underlying DP mechanism differs signif- 



icantly from an ordinary directed percolation process in a porous medium subjected to a 



gravitational field (cf. Sec. 3.9). In ordinary DP the spreading agent represents active sites 



and percolates along the given direction. In the present case, however, water may flow not 
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Fig. 53; Simple model for interface depinning in random media. The pores of the medium are represented 
by cells on a diagonal square lattice. The permeability across the edges of the cells depends on the 
direction of flow: In downward direction all edges are permeable, whereas in upward direction they are 
permeable with probability q and impermeable otherwise. The right panel of the flgure shows a particular 
configuration of open (dashed) and closed (solid) edges. Pumping in water from below, the interface 
becomes pinned along a directed path of solid lines, leading to a finite cluster of wet cells (shaded region). 
The dashed arrow represents an open path in order to illustrate a possible flow. 



only in the direction of the pumping force but also in opposite direction. Moreover, the 
DP process runs perpendicular to the direction of growth, as will be explained below. 



A simple model exhibiting a depinning transition is shown in Fig. 53. In this model 
the pores of the material are represented by cells of a diagonal square lattice. The liquid 
can flow to neighboring cells by crossing the edges of a cell. Depending on the direction of 
flow these edges can either be permeable or impermeable. For simplicity we assume that 
all edges are permeable in downwards direction, whereas in upwards direction they can 
only be crossed with a certain probability q. Thus, by starting with a horizontal row of 



wet cells at the bottom, we obtain a compact cluster of wet cells, as illustrated in Fig. |53. 
The unrestricted flow downwards ensures that the cluster has no overhangs. Clearly, the 
size of the cluster (and therefore the penetration depth of the liquid) depends oxi q. If g is 
large enough, the cluster is infinite, corresponding to a moving interface. If q is sufficiently 
small, the cluster is bound from above, i.e., the interface becomes pinned. 

The depinning transition is related to DP as follows. As can seen in the figure, a 
pinned interface may be interpreted as a directed path along impermeable edges running 
from one boundary of the system to the other. Obviously, the interface becomes pinned 
only if there exists a directed path of impermeable bonds connecting the boundaries of 
the system. Hence the depinning transition is related to an underlying directed bond 
percolation process with probability p = 1 — q running perpendicular to the direction of 



growth. The pinning mechanism is illustrated in Fig. 54, where a supercritical DP cluster 
propagates from left to right. The cluster's backbone, consisting of bonds connecting the 
two boundaries, is indicated by bold dots. The shaded region denotes the resulting cluster 
of wet cells. As can be seen, the interface becomes pinned at the lowest lying branch 
of the DP backbone. Therefore, the roughening exponent coincides with the meandering 
exponent of the backbone 

a = v^v\\ . (256) 

Moreover, by analyzing the dynamics of the moving interface, it can be shown that the 
dynamic critical exponents are given by = a and 5 = 1. Thus, depinning transitions in 
inhomogeneous porous media may serve as possible experimental realizations of the DP 
universality class. 

The prediction ( |256| ) matches surprisingly well with the experimental result a. = 0.63(4) 
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Fig. 54: Pinned interface and the underlying DP process. The figure is explained in the text. 



obtained by Buldyrev et al. p:13| . Therefore, it is near at hand to regard this experiment 
as a first quantitative experimental evidence of DP exponents. However, only one ex- 
ponent has been verified, and it is not fully clear how accurate and reproducible these 
exponents are. Moreover, pressure differences may cause long-range correlations in the 
bulk, leading to a flat interface on large scales. This means that gravity could destroy 
the asymptotic critical behavior. In fact, in subsequent experiments the estimates for the 
critical exponents are scattered over a wide range. For example, Dougherty and Carle 
measured the dynamical avalanche distribution of an air /water interface moving through 



a porous medium made of glass beads p:17|| . According to the DP hypothesis, the distribu- 
tion -P(s) of avalanche sizes s should decrease algebraically. In the experiment, however, a 
stretched exponential behavior P{s) ~ s~^e~^^^ is observed even for small flowing rates. 
The estimates for the exponent b are inconclusive; they depend on the time window of 
the measurement and vary between —0.5 and 0.85. Even more recently Albert et al. pro- 
posed to identify the universality class by measuring the propagation velocity of locally 
tilted parts of the interface |41 . Their results suggest that interfaces propagating in glass 
beads are not described by a DP depinning process, instead they seem to be related to 
the random-field Ising model. Altogether the emerging picture is not yet fully transparent 
and further experimental effort in this direction would be desirable. 

Depinning experiments were also carried out in 2+1 dimensions with a spongy-like 



material used by florists, as well as fine-grained paper rolls |419]. In this case, however, the 
exponent a is not related to 2+1-dimensional DP, instead it corresponds to the dynamic 
exponent of percolating directed interfaces in 2+1 dimensions. In experiments as well as 
in numerical simulations a roughness exponent a = 0.50(5) was obtained. 



6.3 Polynuclear growth 

A completely different type of depinning transition takes place in models for polynuclear 
growth (PNG) [ll4l| , |420H 42l|. A key feature of PNG models is the use of parallel updates, 



leading to a maximal propagation velocity of one monolayer per time step. For a high 
adsorption rate the interface of PNG models is smooth and propagates at maximal velocity 
V = 1. Roughly speaking, this means that the interface is 'pinned' at the light cone of the 
dynamics. Decreasing the adsorption rate below a certain critical threshold, PNG models 
exhibit a roughening transition to a rough phase with v < 1. In contrast to equilibrium 
roughening transitions, which only exist in d > 2 dimensions, PNG models display a 
roughening transition even in one spatial dimension. 
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Fig. 55: Polynuclear growth model. In the first half time step atoms are deposited with probability p. In 
the second half time step islands grow deterministically by one step and coalesce. 



One of the simplest PNG model investigated so far is defined by the following dynamic 



rules 1 141 1- In the first half time step atoms 'nucleate' stochastically at the surface by 



In the second half time step the islands grow deterministically in lateral direction by one 
step. This type of growth may be expressed by the update rule 

hi{t + 1) = max [hi{t + 1/2), hj{t + 1/2)1 , (258) 

ie<i> 

where j runs over the nearest neighbors of site i. 

The relation to DP can be established as follows. Starting from a flat interface hi{0) = 
0, let us interpret sites at maximal height hi{t) = t as active sites of a DP process. 
Obviously, the adsorption process (|257| ) turns active sites into the inactive state with 
probability 1 — p, while the process ( |258D resembles offspring production. Therefore, if 
p is large enough, the interface is smooth and propagates with maximal velocity v = \. 
This situation corresponds to the active phase of DP. Approaching the phase transition, 
we expect the density of sites at maximal height to scale as 



^Y.^h,-t^^P-P-)^^ (259) 



I 

where N denotes the system size and (3 ~ 0.277 is the density exponent of DP. Below a 
critical threshold, however, the density of active sites at the maximum height hi{t) = t 
vanishes after some time, the growth velocity is smaller than 1, and the interface evolves 
into a rough state. Although this mapping to DP is not exact, numerical simulations 
strongly support the validity of Eq. (|25g|). As will be shown below, PNG models are 



actually a realization of unidirectionally coupled DP processes [423, 124 1 



Polynuclear growth models have also been used to describe the growth of colonial or- 



ganisms such as fungi and bacteria [425|. This study was motivated by recent experiments 
with the yeast Pichia membranaefaciens on solidified agarose film |426| . By varying the 
concentration of polluting metabolites, different front morphologies were observed. The 



model proposed in [425| aims to explain these morphological transitions on a qualitative 



level. It is easy to verify that this model follows the same spirit as the PNG model de- 
fined above. They both employ parallel dynamics and exhibit a DP-related roughening 
transition. 

Concerning experimental realizations of PNG models, one major problem - apart from 
quenched disorder - is the use of parallel updates. The type of updates in these models 
is crucial; by using random-sequential updates the transition is lost since in this case 
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Fig. 56: Restricted solid-on-solid growth model exhibiting a roughening transition from a non-moving 
smooth to a moving rough phase. Monomers are randomly deposited whereas desorption takes place only 
at the edges of plateaus. 



there is no maximum velocity. However, in realistic experiments atoms do not move 
synchronously, rather the adsorption events are randomly distributed in time. Therefore, 
random sequential updates might be more appropriate to describe such experiments. It 
thus remains an open question to what extent PNG processes can be realized in nature. 



6.4 Growth with evaporation at the edges of plateaus 

DP-related roughening transitions can also be observed in certain solid-on-solid growth 



processes with random-sequential updates [}427| , |428(| . As a key feature of these models, 
atoms may desorb exclusively at the edges of existing layers, i.e., at sites which have at 
least one neighbor at a lower height. This corresponds to a physical situation where the 
binding energy of atoms at the edges is much smaller than the binding energy of atoms in 
completed layers. By varying the growth rate, such growth processes display a roughening 
transition from a non-moving smooth phase to a moving rough phase. 

A simple solid-on-solid model for this type of growth is defined by the following dynamic 



rules 1 427 1 : For each update a site i is chosen at random and an atom is adsorbed 

hi ^ hi + 1 with probability q (260) 

or desorb ed at the edge of a plateau 

hi mm{hi, hi+i) with probability (1 - q)/2 , 
hi mm{hi,hi-i) with probability (1 — q)/2 . 

Moreover, the growth process is assumed to be restricted, i.e., updates are only carried 
out if the resulting configuration obeys the constraint (|245|) . 



The qualitative behavior of this model is illustrated in Fig. 56. For small q the desorp- 



tion processes ( 261D dominate. If all heights are initially set to the same value, this level 



will remain the bottom layer of the interface. Small islands will grow on top of the bottom 
layer but will be quickly eliminated by desorption at the island edges. Thus, the interface 
is effectively anchored to its bottom layer and a smooth phase is maintained. The growth 
velocity v is therefore zero in the thermodynamic limit. As q is increased, more islands on 
top of the bottom layer are produced until above qc — 0.189, the critical value of q, they 
merge forming new layers at a finite rate, giving rise to a finite growth velocity. 
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Interestingly, this model can be interpreted as a driven diffusion process of two oppo- 
sitely charged types of particles. The charges 

Ci,i+i = h,+l-h^ G {-1,0, +1} (262) 

are bond variables and represent a change of height between two adjacent interface sites. 
In this representation the dynamic rules ( |26[lD -( p6l] ) can be implemented by randomly 
selecting two neighboring bonds and performing the following processes with probabilities 
indicated on the arrows: 

0+ ^ +0 , 00 ^ +- , -0^0-^ _+ ^ 00 . (263) 

{l-g)/2 l-<? {l-g)/2 

In the smooth phase q < Qc, the charges are arranged as closely bound H — dipoles. 
For q > qc, the dipoles become unbound wherefore the fluctuations in the total charge, 
measured over a distance of order A^, diverge with N. Thus the transition can be described 
in terms of correlations between charged particles. 



Relation to directed percolation 

At the transition the dynamics of the model is related to DP as follows. Starting with a flat 
interface at zero height, let us consider all sites with hi = as particles ^ of a DP process. 
Growth according to Eq. (|260| ) corresponds to spontaneous annihilation A — > 0. Con- 
versely, desorption may be regarded as a particle creation process. However, since atoms 
may only desorb at the edges of plateaus, particle creation requires a neighboring active 
site to be present. This process therefore corresponds to offspring production A 2A. 
Clearly, these processes resemble the dynamic rules of a DP process. In contrast to PNG 
models, the DP process takes place at the bottom layer of the interface. Moreover, the 
roughening transition does not depend on the use of either parallel or random-sequential 
updates. 



In the model without the restriction ( p45D , the relation to DP can be proven exactly. In 
this case the processes at the bottom layer decouple from the evolution of the interface at 



higher levels. Introducing local variables Sj = Sh-^ it is possible to map Eqs. ( 260 )-( 261 ) 
exactly onto the dynamic rules 

if Sj = 1 Si ^ with prob. q , 

if {si-i,Si, Si+i} = {0, 0, 1} Si ^ 1 with prob. (1 - q)/2 , 

if {si_i, Sj, Sj+i} = {1, 0, 0} Si^ 1 with prob. (1 - q)/2 , 

if {si_i,Si,Si+i} = {1,0, 1} Si^l with prob. 1 - g . 



These rules define a contact process on a square lattice (cf. Fig. in which particles 
self-annihilate at unit rate and create offspring with rate A/2 = (1 — q)/2q. Since the 
percolation threshold of the contact process Ac = 3.29785(8) is known, we can predict the 
critical point of the unrestricted growth to be given by qc = 0.232675(5). model. For the 
restricted model there is no exact mapping to a contact process. Nevertheless numerical 
simulations strongly suggest that the bottom layer still exhibits DP behavior. 

The underlying DP transition is the origin of the roughening transition. In the active 
phase of DP the interface fluctuates close to the bottom layer so that the interface is 
smooth. Approaching criticality the bottom layer occupation no vanishes as 

no ~ {qc - qf , (264) 
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Fig. 57: Spontaneous symmetry breaking. Left: Simulation in the smooth phase starting from a random 
interface configuration. Black (white) pixels indicate even (odd) heights. As time evolves the systems 
'coarsens', providing a robust mechanism for the elimination of minority islands. Right: Order parameters 
M„ versus time at criticality starting from a flat interface. 



where /3 is the density exponent of DP. In the inactive phase of DP the interface detaches 
from the bottom layer and evolves into a rough state. Therefore, the front velocity v for 
g > (7c is proportional to the characteristic survival time of a DP process in the inactive 
phase 

V ~ l/^ii ~ (g - q.fl-W . (265) 

With respect to the microscopic rules the roughening transition in these models seems to 
be as robust as a DP transition. Including the dynamics at higher levels of the interface, 
the models turn out to be described by unidirectionally coupled DP processes (see below). 



Scaling of the interface width 

Turning to the scaling properties of the interface width (246) the emerging picture is less 
clear and may indicate the presence of many length scales. With only a single length scale 
one would expect the saturation width in finite systems to scale as Wgati]^) ~ in the 
growing phase and Wsati^) ~ at criticality, where a and a' are generally different 
critical exponents. Thus the expected scaling form would read 



Wsat{N,q)=N^' g{N'/''Hq-qc)) , 
where g{u) is a universal scaling function with the asymptotic behavior 



\u\ " ' 
const 



\u 



{a—a')u± 



for u — oo 
for u ±0 , 
for u +00 



(266) 



(267) 



This scaling form implies that the width in a finite size system saturates at 

'(g^_^)-"W if q<qc, 

Wsat{N,q) = Im' ifq = qc, 

N^iQc-qY"-"'^"^ iiq<qc- 



(268) 
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However, the width at criticahty actually increases as w{t) ~ (Inty until it saturates at 

^l;,,^~(lniV)^ (269) 

where k ~ 0.43, suggesting that a' = 0. Moreover, as the growth rate is increased 
and the interface starts to move, the width diverges as Wsatig) — WsatiQc) ~ — Qc)^'^^, 
suggesting that a ~ 0.9. This value differs from the expected value a = 1/2 for growing 
interfaces in one spatial dimension. Therefore, the scaling behavior of the interface width is 
presumably characterized by many different length scales. This point of view is confirmed 
by a numerical analysis of correlation functions at different levels |428]. 



Spontaneous symmetry breaking 

The roughening transition in growth models with evaporation at edges of islands is accom- 
panied by spontaneous symmetry breaking of translational invariance in height direction. 
This symmetry is associated with a family of non-conserved magnetization-like order pa- 
rameters 

The order parameter Mi = | 1)'*^ |, for example, measures the difference between 
the densities of sites at even and odd heights, respectively (see Fig. In the smooth 
phase the order parameters M„ tend towards stationary positive values. Approaching the 
phase transition these values vanish algebraically as 

(M„) ~ {q, - qf- , (271) 

where 6i ~ 0.65, 02 ~ 0.40, and ^3 ~ 0.23. At criticality, starting from a flat interface, the 
order parameters decay as Mn{t) ~ t~^"^'^^^, as shown in the right hand graph of Fig. 57. 



In the rough phase q > qc, where the heights hi become uncorrelated over large distances, 
Mn = in the thermodynamic limit. 

Experimental realizations 

With respect to experimental realizations of the dynamic rules ( |260| )-( |26l| ) it is important 
to note that atoms are not allowed to diffuse on the surface. This assumption is rather 
unnatural since in most experiments the rate for surface diffusion is much higher than 
the rate for desorption back into the gas phase. Therefore, it will be difficult to realize 
this type of homoepitaxial growth experimentally. However, in a different setup, the 
above model could well be relevant p:29(| . As illustrated in Fig. ^ a laterally growing 
monolayer could resemble the dynamic rules ( |260D and ( |261| ) by identifying the edge of the 
monolayer with the interface contour of the growth model. In this case 'surface diffusion', 
i.e. diffusion of atoms along the edge of the monolayer, is highly suppressed. Moreover, 
in single-step systems (such as fcc(lOO) surfaces) it would also be possible to implement 



the restriction (245) 



6.5 Unidirectionally coupled directed percolation processes 

So far we have seen that the scaling properties of quantities involving only the bottom 
level of the interface may be adequately described by using DP exponents. In particular it 



121 



Fig. 58: Possible experimental realization of the growth process defined in Eqs. (260) and (261) by lateral 
growth of a monolayer on a substrate. 




Fig. 59: Coupled directed percolation. Stationary densities Uk as functions of the distance from criticality 
in (a) the unrestricted growth model, (b) the restricted growth model, and (c) a unidirectionally coupled 
sequence of directed bond percolation processes. 

was shown that the density of exposed sites at the bottom layer decreases as uq ~ {qc — q)^ 
until it vanishes at the transition. Let us now turn to the scaling properties of the first 
few layers k = 1,2, . . . above the bottom layer. Since the scaling properties at the bottom 
layer k = m the smooth phase are completely characterized by the three DP exponents 
/3, z^_L, and it is natural to assume that the next layers obey similar scaling laws with 

analogous exponents, /J'^^-', and t^^l'^ ■ For example, the densities 

1 ^ 

"*^ = iV^^'^^-^' ' k = 0,l,2,... (272) 

j=0 i 

of sites at height < k above the bottom layer are expected to scale as 

nk--{qc-qf'^ ■ k= 1,2,3,... (273) 

In principle all these exponents could be different and independent from each other. How- 
ever, extensive numerical simulations and field-theoretic considerations [ |42^ ] suggest that 

(k) (k) 

the scaling exponents and are identical on all levels and equal to the DP ex- 
ponents z^_L and This remarkable property implies that the growth process at criti- 
cality is characterized by a single dynamic exponent z = The density exponents 
^(1) ^ ^(2) ^ ^(3) ^ , _ ^ however, turn out to be different and considerably reduced compared 
^(0) = p_ These exponents appear to be non-trivial in the sense that they are not simply 
related to DP exponents. 

In order to explain the reduced values of P^''^ it is useful to introduce the following 
particle interpretation. Assuming that the minimal height of the interface is zero, lattice 
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sites at height /ij < 0, 1, 2, . . . may be associated with particles A, B,C, . . . , respectively. 
By definition, particles of different species are allowed to occupy the same site simulta- 
neously. For example, the presence of an A-particle induces simultaneous presence of all 
other particle species at the same site. As shown before, the A-particles of the unrestricted 
model evolve independently according to the dynamic rules of a contact process. Similarly, 
in absence of A-particles, the -B-particles will evolve according to the rules of a contact 
process. In presence of an ^-particle, however, 5-particles are instantaneously created 
at the same site, giving rise to an effective reaction A ^ A + B at infinite rate. As this 
reaction does not modify the configuration of the A-particles, it couples the two subsys- 
tems in one direction without feedback. Similarly, the -B-particles induce the creation of 
C-particles by an effective reaction B B + C. Thus we obtain a unidirectionally coupled 
sequence of contact processes corresponding to the reaction-diffusion scheme 

A^2A B ^2B C ^2C 

A^ A + B B ^ B + C C ^C + D etc. 

Notice that this sequence can be truncated at any level without changing the dynamics of 
the lower levels. In fact, numerical simulations of the growth model and a unidirectionally 
coupled sequence of three (l+l)-dimensional directed bond percolation processes yield 
compatible estimates for the critical exponents (cf. Fig. 

/?(°) = 0.28(1) , = 0.13(2) , = o.05(l) . (274) 

Thus, 'coupled DP' is a universal phenomenon and should play a role even in more gen- 
eral contexts, namely whenever DP-like processes are coupled unidirectionally without 
feedback. 

Mean field approximation 

The simplest set of Langevin equations for unidirectionally coupled DP reads |l428|| 

at0fc(x, t) = KkM^, t) - A(/.i(x, t) + DV^^ki^, t) + ^(/)fc_i(x, t) + Cfc(x, t) , (275) 
where Ck are independent multiplicative noise fields with correlations 

(a.(x,t)) = 0, (Cfc(x,t)C/(x',0) = 2r0fc(x,t)5fc,iJ'^(x-x')5(t-O • (276) 
Assuming that 0_i = 0, the lowest equation for k = reduces to the ordinary Langevin 



equation of DP ( 139 ). The parameters Kk control the rates for offspring production at 
level k. In principle all these parameters could be different, leading to interesting mul- 
ticritical behavior [423]. Here we will restrict to the simplest case where the parameters 



coincide. 

The reduced values of P^^^ , /J'-^-* , . . . can already be understood on the level of a sim- 
ple mean field approximation. Determining the stationary solutions of the mean field 
equations 

dt(l)k = n4>k - ^4>k + /^0fc-i (277) 

and expanding the result for small values of k, the fields (pt are found to scale asymptoti- 
cally with the mean field exponents 

(5^Sf = 2"' • (278) 
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In addition, dimensional analysis of Eqs. ( ^75| )-( p76| ) reveals that the mean field scaling 
exponents u^'^ and i^y*^^ coincide with the DP exponents = 1/2 and = 1 on all 
levels, i.e., they do not depend on k. 

Field-theoretic treatment 

The above mean field approximation is expected to hold above the critical dimension 
dc = 4. In less than four dimensions fluctuation corrections to f3^^^ have to be taken into 
account. These corrections can be approximated by field-theoretic renormalization group 



techniques 1 425, 424 1. To this end the master equation is mapped onto a second-quantized 



bosonic operator representation, leading to the effective action (cf. Sect 

K-l 



S[4>o,(t>Q,4>i,4>i,- ■ ■]= / d'^xdt l4)i,{x,t) (rdt- DV"^ - Kj (j)k{x,t) - fi(f>k(j)k-i 

k=o '' I 

+ ^'^/c(x,t)(^</)fc(x,t) - (^fc(x,t)^(/'fc(x,t)| , 

(279) 

where K is the total number of levels in the hierarchy. Because of the unidirectional 
structure the truncation of the hierarchy at finite K does not affect the temporal evolution 
of the subsystems k < K. The field-theoretic treatment turns out to be rather complex, 
even a one-loop calculation for a two-level system involves 34 different diagrams. Details 



of these calculations are given in Ref. |424], whose main results we summarize here. A 



one-loop calculation in the inactive phase reveals that the RG flow is characterized by 
two flxed lines. One of them is unstable and corresponds to a situation where the two 
systems are decoupled. The other one is stable and describes the interacting case. One 
can show that strongly ultraviolet-divergent contributions cancel along the stable fixed 
line. To determine the exponent similar calculations have to be performed in the 
active phase. Choosing a particular point along the fixed line it is possible to derive the 
result 

/?(0)=i_e/6, /3W = l/2-e/8, (280) 

where e = 4 — d. This result explains the downward correction of the critical exponent 
P^^^ in less than four dimensions. 

The field-theoretic analysis involves various technical and conceptual problems. On the 
one hand, in the active phase several infrared-divergent diagrams are encountered [}43Cl(| , 
without being clear to what extent they will affect the physical properties of coupled DP. 
On the other hand, the coupling constant n between different levels is shown to be a 
relevant quantity, i.e., it grows and finally diverges under RG transformations. This may 
be the cause for numerically observed violations of scaling. In fact, the curves in Fig. ^ 
for k > 1 are not perfectly straight but slightly bent. This curvature is neither related 
to transients nor to finite size effects. A careful analysis shows that the field-theoretic 
prediction seems to apply in an intermediate scaling regime, whose size depends on the 
coupling constant fi. The breakdown of scaling may be an artifact of the lattice realization 
which, in contrast to the field-theoretic prediction, seems to limit the value of fj,. 
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Fig. 60: Roughening transition in a model for dimer adsorption and desorption. Left: Dynamic rules for 
adsorption and desorption at the edges of terraces. Right: Typical interface configurations below, at, and 
above criticality. 



6.6 Parity-conserving roughening transitions 

So far we discussed two examples of depinning transitions related to DP, namely interface 
depinning in driven random media, polynuclear growth, and growth without evaporation 
at the edges of terraces. It is therefore interesting to ask whether it is possible to find 
examples of roughening transitions with non-DP behavior in one spatial dimension. 



As discussed in Sec. 4.2 non-DP phase transitions into absorbing states can be observed 
in systems with additional symmetries. For example, non-DP behavior is observed in 
parity- conserving branching processes. Therefore, the question arises whether it is possible 
to replace the underlying DP transition in the previously discussed growth models by a 



parity-conserving mechanism. As shown in Ref. |431], this can be done by considering 
the growth of dinners with evaporation at edges of plateaus. As dimers consist of two 
atoms, the number of particles at each height level is preserved modulo 2. This definition 
of the dynamic rules mimics a BAWE at the bottom layer of the interface. It turns out 
that some of the critical properties of the roughening transition in this model are indeed 
characterized by PC exponents (cf. Eq. ( ^2l[) ). 



As shown in Fig. the model generalizes the growth model introduced in Sec. 6.4, 
replacing monomers by dimers. In each attempted update two adjacent sites i and i + 1 
are selected randomly. If the heights hi and /ij+i are equal, one of the following moves is 
carried out. Either a dimer is adsorbed with probability p 

hi^hi + 1, hi+i hi+i + 1 (281) 
or desorbed with probability 1 — p: 

hi, hi+i min(/ij_i, hi, /ij+i, /ij+2) • (282) 
An update will be rejected if it leads to a violation of the RSOS restriction ( ^45| ) . 



The transition of the model is illustrated in Fig. 60: If g is very small, only a few 
dimers are adsorbed for a short time so that the interface is smooth and pinned at the 
(spontaneously selected) bottom layer. As q increases, a growing number of dimers covers 
the surface and large islands of several layers stacked on top of each other are formed. 
When q exceeds the critical value qc — 0.317(1), the mean size of the islands diverges and 
the interface evolves into a rough state. 
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Fig. 61: Dimer model: (a) Finite size scaling of the width, (b) Decay of the densities no,ni,n2 as a 
function of time, (c) Decay of the corresponding densities in unidirectionally coupled BAWE's. 



Naively one may expect the interface to detach from the bottom layer in the rough 
phase, resulting in a finite propagation velocity. In the present case it turns out, how- 
ever, that the interface remains pinned to the initial height, i.e., it does not propagate 
at constant velocity. This is due to the fact that a stochastic deposition process cannot 
create a dense packing of dimers. Instead the emerging configurations are characterized 
by a certain density of defects (solitary sites at the bottom layer) where dimers cannot 
be adsorbed. Because of the RSOS condition ( |245| ) these defects act as 'pinning centers' 
which prevent the interface from growing, leading to triangular configurations as shown 



in Fig. |60|. The pinning centers cannot disappear spontaneously, they can only diffuse by 
interface fiuctuations and recombine in pairs so that their number is expected to decrease 
very slowly. Therefore, the interface of an infinite system in the rough phase does not 
propagate at constant velocity. Instead the squared width and the average height grow 
logarithmically with time as h{t) ~ w^t) ~ logt. At criticality, this behavior crosses over 
to w'^{t) ^ h{t) ~ \ogt. Thus the expected scaling form reads 

t) ~ log \t-'^ F{t/L^) \ , (283) 

where -F is a universal scaling function and 7 plays the role of a growth exponent. For 
7 = 0.172(10) and z = 1.76(5) we obtain a fairly accurate data collapse (see Fig. |6l|a), 
supporting the supposition that z is the dynamic exponent of the PC universality class. 
Similarly the densities (see Eq. ( ^721) ) are found to decay at criticality as ni^{t) ~ t~^^ 
with 



(5o = 0.280(10) , 5i = 0.200(15) , 82 = 0.120(15) , (284) 

where 5q = f3/v\\ is the usual cluster survival exponent of the PC class. 

Using the particle interpretation of Sec. the temporal evolution of the A-particles 
resembles a BAWE. Similarly, the S-particles perform an effective BAWE on top of inac- 
tive islands of the A-system. Since A-particles instantaneously create i3-particles, the two 
subsystems are coupled by the effective reaction A ^ A + B at infinite rate. As this reac- 
tion does not modify the configuration of the A-particles, it couples the two subsystems 
only in one direction without feedback. On the other hand, the RSOS condition ( p45D 
introduces an effective feedback so that the A-particles are not completely decoupled from 
the 5-particles. However, it seems that the inhibiting influence of the S-particles does not 
affect the critical behavior of the 74-particles. Similarly, the C-particles are coupled to the 
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Fig. 62: Wetting of a surface - schematic illustration. 



i?-particles by the effective reaction B ^ B + C . Therefore, the dimer model resembles a 
unidirectionally coupled sequence of BAWE's, corresponding to the reaction scheme 

A^2,A B ^3B C ^3C 

2A^0 2B^0 2C ^ 

A^ A + B B ^ B + C C ^C + D etc. 



In fact, as shown in Fig. 61 b-c, a coupled hierarchy of three BAWE's shows the same 
critical behavior as the dimer model at the first few layers. 



6.7 Nonequilibrium wetting transitions 

Wetting phenomena are observed in various physical systems where a bulk phase (e.g. a 
gas or liquid) is brought into contact with a wall or a substrate. Because of interactions 
between bulk phase and surface, a thin layer of another phase may be formed which is 
attracted to the substrate. The thickness of the layer fluctuates and may depend on various 
parameters such as temperature or chemical potential. As some of these parameters are 
varied, the thickness of the layer may diverge, leading to a wetting transition. Theoretical 
models for wetting phenomena ignore the details of molecular interactions between surface, 
layer and bulk phase. Instead, they characterize the system by an interface without 
overhangs that separates the two phases. The configuration of the interface is given in 
terms of the height h{x) > of the interface at point x on the surface. Within this 
approach, wetting transitions may be viewed as the unbinding of the interface from the 
wall. 

Wetting transitions at thermal equilibrium have been theoretically studied and exper- 
imentally observed in a large variety of systems (for a review, see Ref. [ |432| ). Models for 
equilibrium wetting are usually defined by an effective Hamiltonian of the form 

n = J d'^-'x [^(Vhf + V[hix)]] , (285) 

where a denotes the surface tension of the interface and d — 1 the dimension of the in- 
terface [ [433| - |435| . The potential y[/i(x)] yields the effective interaction between wall and 
interface. Usually it contains an attractive component binding the interface to the wall. 
However, as temperature or other parameters are varied, the attractive component of the 
potential may become so weak that the potential is longer able to bind the interface. 
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Fig. 63: Phase diagram of the nonequilibrium wetting model. The wetting transition takes place along 
the solid line. Along the dashed line the growth velocity in the model without a wall does not depend on 
a global tilt of the interface, indicating that the effective KPZ nonlinearity vanishes. Both lines intersect 
in the point q — p = 1, where the microscopic processes obey detailed balance. 



leading to a wetting transition. In d = 2 dimensions one usually distinguishes between 
critical and complete wetting. Critical wetting refers to the divergence of the interface 
width when the wetting transition is approached by moving along the coexistence curve of 
the bulk and surface phases. On the other hand, complete wetting refers to the divergence 
of the interface width when the chemical potential difference between the two phases is 
varied, moving towards the coexistence curve. Critical and complete wetting transitions 
are associated with generally different critical exponents. 

While equilibrium wetting transitions are well understood, the investigation of wetting 
transitions under nonequilibrium conditions has started only recently. Here the surface 
layer is adsorbed to the wall by a growth process whose dynamics, unlike in equilibrium 
processes, does not obey detailed balance. As expected, the resulting critical exponents 
differ from those observed at thermal equilibrium. 



A lattice model for nonequilibrium wetting 

A simple lattice model for nonequilibrium wetting can be defined as follows [436|. As in 



Sec. 6.4, the interface is given by height variables /ij = 0, 1, . . . , oo. For each update a site 



i is selected randomly and one of the following moves is carried out: 

/ij — > /ij + 1 with prob. q/{q + p + I) , 

hi min(/ij_i, hi, /ij+i) with prob. l/{q + p + 1) , (286) 

hi ^ hi - 1 if = hi = /ij+i with prob. p/{q + p + 1) ■ 

The selected move will be rejected if it would result in a violation of the RSOS constraint 
\hi — /ij+il < 1. In addition, a hard-core wall at zero height /i = is introduced, i.e., a 
process is only carried out if the resulting interface heights are non-negative. Generally 
the processes defined above do not satisfy detailed balance. 
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By varying the relative rates of these processes, a transition from a binding to a non- 
binding phase is found. This wetting transition can be understood as follows. Without the 
wall for fixed p > 0, the parameter q controls the mean growth velocity of the interface. 
This velocity may be positive or negative and vanishes at some critical value q = qc- On 
large time scales a lower wall will only affect the interface dynamics if the interface moves 
downwards, i.e. q < qc, leading to a smooth interface. In the moving phase q > q^ 
however, the interface does not feel the wall and evolves into a rough state (see Fig. |63| ). 
It is obvious that in the moving phase the interface velocity scales as w ~ (g — qc)^ with 
y = 1. In the smooth phase q < qc, the expected scaling for bottom layer occupation and 
width is given by 

PO ~ {qc - qT" , w^iqc- q)'^ , (287) 

where xq and 7 are certain critical exponents. The above model includes two special 
cases, namely p = 0, where the interface cannot move below its actual minimum height, 
and p = 1, where the dynamic rules satisfy detailed balance. 



The exactly solvable case p = 1 

For p = 1 the dynamic rules ( |286| ) can be mapped onto an exactly solvable equilibrium 
model which exhibits a transition to complete wetting. For q < 1 the probability of finding 
the interface in the configuration {hi, . . . , /iat} is given by the distribution 

P{hi,... ,hN) = P{H) = Z-^ q^^^^'-'^^\ (288) 

where H = H{hi, . . . , h^) = '}2d=i heights and Z = J2hi hjv 

denotes the partition sum running over all possible interface configurations. Eq. (|28^ ) can 
be proven by verifying the detailed-balance condition 

w{aH (JH+i)/w{aH+i (Th) = q. (p = 1) (289) 

which is consistent with the hard wall constraint /ij > 0. The steady state distribu- 
tion ( ^881) does not exist for q > 1 where the interface propagates at constant velocity. 
The critical exponents xq and 7 can be computed by analyzing the transfer matrix acting 
in spatial direction [437, 43^ 



r,,, = l^' if|/i-/i'|<i , 

\ otherwise ' 

where h, h' > 0. Steady state properties can be derived from the eigenvector (p that 
corresponds to the largest eigenvalue /i of the transfer matrix 'Yl'h'=o'^h,h'4>h' = P'4>h- 
From the squares of the eigenvector components various steady state quantities can be 
derived. For example, the probability ph of finding the interface at height h is given by 
Ph = 'Ph/ Ylh' ^h'- Close to criticality, where e = 1 — q is small, one can carry out the 
continuum limit (j)h — > (p{h), replacing the discrete heights h by real- valued heights h. 



Then, the above eigenvalue problem turns into a differential equation [437| which, to 
leading order in e, is given by 



+ {3- n)- Sehj (p{h) = . (291) 
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Simple dimensional analysis indicates that the height variables scale as /i ~ g~i/3 wherefore 
the width diverges as uP ~ e~^/^. Similarly one can show by elementary calculations that 
po ~ £• The critical exponents for p = 1 are thus given by 

xo = l, 7 = 1/3. (292) 



The generic case p ^ 1 

In the nonequilibrium case p ^ 1 the critical exponents can only be approximated numer- 
ically. For p = 0, where the interface cannot move below its actual minimum height, the 
hard-core wall becomes irrelevant and the model reduces to a growth process similar to 



the one discussed in Sec. S^, belonging to the universality class of unidirectionally coupled 
directed percolation. For p = 1 the numerical estimates are consistent with the previously 
derived exact results. For < p < 1 a very slow crossover to a different behavior with 
critical exponents y = 1.00(3), xq = 1.5(1), and 7 = 0.41(3) is observed. Similar results 
are obtained for p > 1 except for xq which is close to 1 in this case. 

It is believed that this type of nonequilibrium wetting can be modeled by the KPZ 
equation (|252| ) with an additional term for the effective interaction between the wall and 
the interface ||438|,E4C 



= „ + .V'Hr, t) + A(VMr. *))^ + C(r, t) - . (293) 

ot on{r, t) 

Dimensional analysis suggests that the width exponent described by this equation should 
be given by 7 = (2 — z)/{2z — 2), where z = 3/2 is the dynamic exponent of the KPZ 
universality class, yielding 7 = 1/2. The numerical estimates of 7 ~ 0.41 are smaller, 
presumably because of the very slow crossover to the exactly solvable case 7 = 1/3. In 
addition, Eq. ( |293| ) suggests that the bottom layer occupation po should be proportional 
to the inverse correlation length, hence xq = i'± = l/(2z — 2) = 1. However, this scaling 
argument seems to hold only for p > 1, whereas for < p < 1 much larger values for xq are 



observed. As shown in Ref. [436|, the changing sign of the coefficient A in Eq. ( ^931) leads 
to different universality classes in both cases, corresponding to the distinction between an 
'upper' and a 'lower' wall in Ref. [}440|| . In fact, comparing the growth velocities of a flat 
and a tilted interface in absence of a wall it is found that the nonlinear term (V/i(r, t))^ 
is indeed non- vanishing in the {p, q) plane, except for the dashed line shown in Fig. |6^. 
Therefore, moving along the transition line, the sign of A changes precisely at the integrable 
point, leading to a different exponent xq. 

Nonequilibrium wetting of an attractive substrate 

In the above wetting model the substrate is introduced as a hard wall at zero height. 
Therefore, the model neglects interactions between the substrate in the surface layer. 
Loosely speaking, the free energies of the exposed and the wetted substrate are assumed 
to be equal. In order to describe more realistic situations, the model has to be generalized 
by taking interactions between the substrate and the surface layer into account. This 



can be done by introducing a modified growth rate qq at zero height |441]. Obviously the 
attractive short-range interaction at the bottom layer is a surface effect. Therefore, the 
critical point Qc remains unchanged. However, if the interaction is strong enough, the 
transition becomes discontinuous. In the equilibrium case p = 1 it can be proven by using 
transfer matrix methods that for go < 2/3 there is a first-order wetting transition. 
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Fig. 64: Phase coexistence in a nonequilibrium wetting model. Left: Phase diagram of the wetting model 
with an attractive substrate qo — 0.2. Right: Elimination of an island in the coexistence phase. The 
flat island (a) grows quickly until it reaches a triangular shape (b). Because of the KPZ nonlinearity, 
the deposition rate decreases, leading to a laterally shrinking island (c). Finally the island is eliminated, 
ensuring phase coexistence. 



In the nonequilibrium case the morphology of the phase transition depends on the 
sign of the KPZ nonlinearity. For p > 1 the emerging picture is essentially the same 
as for p = 1, although with different critical exponents. For p < 1, however, there is a 
whole region in the phase diagram where the pinned and the moving phase coexist. As 
illustrated in Fig. islands generated by fluctuations quickly grow until they reach an 
almost triangular shape. Since the KPZ nonlinearity is negative, adsorption processes at 
the inclined edges of the islands are strongly suppressed, allowing the attractive interaction 
to reduce the size of the island in lateral direction. This ensures the stability of the pinned 
phase in parts of the phase diagram where a free interface would move away from the wall. 
This model demonstrates that nonequilibrium effects do not only lead to different critical 
exponents but may also change the whole phase structure of a model. 
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Appendix A: Vector space notation and tensor products 



A one-dimensional lattice model, whose sites i = 1, . . . ,N are either occupied (sj = 1) 
or vacant (sj = 0), can be in 2^ different states s = {si, S2, ■ ■ ■ , sn}. In order to represent 
the probability distribution Pt{s) as a vector in a 2^-dimensional vector space let us define 
an orthogonal set of basis vectors \s) corresponding to the configurations of the system. 
Using the representation 



|0) 



ID 



the basis vectors are given by 

\s) = \si) (g) |S2) (8) ... (8) \sn) , 

where '(8)' denotes the tensor product of two vectors: 



(A.l) 



(A.2) 



/'aibi\ 
\a2b2J 



(A.3) 



Row vectors {s\ are the transposed vectors of |s). In this vector space the probability 
distribution Pt{s) can be represented by the vector 



s s 



(A.4) 



Defining the sum vector over all states 



(1,1,... ,1) 



(A.5) 



the normalization of the probability distribution can be simply expressed as {l\Pt) = 1- 
Similarly the ensemble average {A{t)) of any observable A can be expressed as 



m)) = {i\A\Pt] 



(A.6) 



The empty lattice is represented by the vector \vac) = |0)®^. Local operators act only on 
a finite number of adjacent sites. For example, a single-site operator Ai can be written as 

= 1 (g) 1 . . . ^A^ . . . (g) 1 , (A. 7) 

i-th position 

where 1 and A are 2x2 matrices. Here the tensor product of two matrices is defined by 



ai a2 
03 a4 







bi 
63 



b2 
64 



f aibi 0162 02^1 a2&2\ 

0163 0164 0263 02^4 

0361 0362 04^1 04^2 

\a363 0364 0463 0464/ 



(A.^ 



Similarly one can define two-site operators where is a 4 x 4 matrix. The above 

notations can be easily generalized to systems with n > 1 particle species by introducing 
local vectors with n components in Eq. ( |A.l| ). 
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Appendix B: Derivation of the effective action 



In order to derive the effective action of Reggeon field tlieory by integration of the 
noise, we first introduce a response field (/'(x, t). This allows the (5- function to be expressed 
as an oscillating integral 



Z 



DC P[C] / D4) 0] exp \i / d'^x dt 4>{dt(l) - DV^cj) - K(p + Xcf? - C) 



(B.l) 



where (j)] denotes the Jacobian. After a Wick rotation in the complex plane the 
dependent noise contribution can be separated by 



D(l) 4>] exp - J d'^x dt ^{dtcj) - OV'^cp - K<p + Xcp'^) 

DCP[(] expf / d'^xdt 0C 



Z 



Because of the correlations (|140| ) the probability distribution P[C] is given by 

C'(x,t) 



(B.2) 



P[C] = m exp 



d'^xdt 



(B.3) 



2r</.(x,t), 

where f[(j)] is a (field-dependent) normalization factor. This allows the noise to be inte- 
grated 



DCP[C] exp / d'^xdt 0C 



DC, exp 



DC exp 



d'^xdti^C 



2T<i) 



d'^xdti-^^^ 



2T(t) 



(B.4) 



m exp(^ / d^xdt -^^2 



where we used Gaussian integration of the form 

1 . .2 



+ 00 

drj 

-oo v27rC(/> 



exp 



(j)ri = expf- 



The resulting partition function reads 



Z D(t>D4>I'\ 



exp 



d'^x dt (j) [dtc/} - DV^ 



V2 



k6 + 



(B.5) 



(B.6) 



It is convenient to symmetrize the cubic terms in the partition function. To this end we 
rescale the fields by 



V = V2fX 



(B.7) 



In order to keep the action symmetrized during the RG procedure, one has to introduce 
an additional coefficient r in front of the time derivative. Dropping the primes the effec- 
tive action S = Sq + Sint is given by the expressions ( |144| ) and ( |145| ). Alternatively the 
action may directly be derived from the master equation of a contact process by introduc- 
ing bosonic creation and annihilation operators (for this standard procedure we refer to 
Refs. m, 11,111). 
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d 


1 


2 


3 


4 


5 


6 




2 


27r 




27r^ 


87rV3 




Kd 


j_ 

7r 


1 

27r 


1 


1 

8^ 


1 


1 

647r^ 


Vd 


2 


vr 


47r/3 


^V2 







Table 5: Surface area Sd, Kd = Sd/i^iv)'', and the volume Vd of a d-dimensional sphere. 



Appendix C: Shell integration 



The one-loop integrals in Wilson's renormalization group approach take the form 



1 



{2ttY 



<fk' f{k',hk',k'^) 



(C.l) 



where '>' denotes integration in the momentum shell $7(1 — I) < \k'\ < Q . This integral 
can be written as 



(2^ 



Sd-i [ de sm'^-^ 9 f{k'^, n \k\ cos 6, Q^) , (C.2) 

^0 



where Sd and Vd denote surface area and volume of a d-dimensional sphere: 



S, 



2tt^ 
T{d/2) 



T/ _ 



(C.3) 



We also use the notation Kd = Sd/i^n)'^. For easy reference we listed some of the values 
for Sdi Kd, and Vd in Table ^. 

To evaluate Eq. (C^) it is often helpful to use the formulas 

r d9 sin''-' 9 = , 

Jo '-'(i-l 

r d9 sin'^-^^cos^^ = -1^ . 
Jo dbd-i 

In particular, if the function / does not depend on 9 the integral I{k) reduces to 



(C.4) 
(C.5) 



I{k) = £ d'^k' f{k\ k'^) = IKd^'' fik^n'). (C.6) 
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Appendix D: One- loop integrals for directed percolation 

The integral for propagator renormalization reads 

= I Dk'u'Go{'l + k',^+u;')Go{'^-k',^-u') (D.l) 



> 

1 



Dk'io 



> + -K- iT{^ + u;')) (^(1 - k'y -k- iT{^ - u')') 

Denoting A"^ = ± k')"^ — k — ir^ and integrating with respect to the pole irco' = A'^ 
we obtain 

tP = ^ / rldj,' ^' \ 

(27r)'^+i y> ' ^ {A+ - IT J) (A" + IT J) 



(27r)'^+iir y> + A- 

2{27ryTjy \Dk'^ + Dk'^ 



2T(\Dk'^ + n^D -K- iro;) 
The integral for vertex renormalization is given by 
= J DkuGl{k,u)Go{-k,-u) 



f oo 



1 ' jdj 



(fk / duj-—^ : : . (D.3) 



(27r)<^+i 7> [Dk'^ -K- iTuf {Dk'^ -n + iru) 

Integration with respect to the poles iru = Dk^ — k yields the result 
rV 27ri r 1 iK^n^ 



= — / d'^k = = ° fD 4) 
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Appendix E: Notations 



Frequently used symbols 



A,B,C,... 


particle species 





vacant site 


a 


roughness exponent 


P 


density exponent 




density exponents at different hierarchy levels 


S 


exponent for temporal decay 


D 


diffusion constant 


r 


noise amplitude 


d 


spatial dimension of the system 


D,E,D,E 


matrices for matrix product states 


A 


distance from criticality 




Hamming distance (damage) 


7 


growth exponent 


r 


noise amplitude 


hi{t) 


interface height at site i 


n 


Hamiltonian, internal energy 


< i > 


set of nearest neighbors of site i 


e 


critical initial slip exponent 


On 


critical exponents of the order parameters M„ 


h 


probability for empty interval of length £ 


J 


coupling constant in equilibrium systems 


K 


rate for offspring production in Langevin equations 


L 


length of the system 


X 


coefficient for nonlinear term in Langevin equations 


p(x, t) 


particle density 


A 


dilatation parameter for scaling transformations 




coupling constant between different levels 


C 


Liouville operator 


Mn 


order parameters for spontaneous symmetry breaking 


nk{t) 


density of sites at height k above the bottom layer 


N 


system size (total number of sites N = L'^) 




temporal scaling exponent 




spatial scaling exponent 


P 


percolation probability 




probability to find the system in state s at time t 


m) 


local persistence probability 


Pgit) 


global persistence probability 


q 


probability for interface growth 


s 


state of the model 


Si 


local state Sj = 0, 1, . . . of lattice site i 
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s 


field-theoretic action 


a 


control exponent for anomalous DP 




local Ising spin cTj = ±1 at site i 


T 


temperature 


V 


interface velocity 




transition rate from state s to state s' 




interface width 


n 


cutoff in momentum space 


C(x,t) 


noise field in Langevin equations 




spatial correlation length 


^11 


temporal correlation length 


z 


dynamic critical exponent 


Z 


partition sum 



Abbreviations 



BAWE 


branching-annihilating random walk with even number of offsprinj 


CDP 


compact directed percolation 


DK 


Domany-Kinzcl (model) 


DMRG 


density matrix renormalization group 


DP 


directed percolation 


DP2 


directed percolation with two absorbing states 


DS 


damage spreading 


EW 


Edwards- Wilkinson 


IMF 


improved mean field 


IPDF 


interparticle distribution function 


KPZ 


Kardar-Parisi-Zhang 


MC 


Monte Carlo 


MF 


mean field 


MPS 


matrix product state 


PNG 


polynuclear growth 


RG 


renormalization group 


RSOS 


restricted solid on solid 


SOS 


solid on solid 
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